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Chapter  1 

Introduction  and  Review 
of  Current  Methods 

1.1  Background 

Wave  kinematics  are  involved  in  most  processes  that  coastal  engineers  study. 
Knowledge  of  fluid  velocities  and  accelerations  are  necessary  for  the  study  of  the 
wave  loading  of  structures  through  the  use  of  the  O’Brien-Morison  equation.  Knowl¬ 
edge  of  the  kinematics  near  the  sea  bed  are  necessary  for  studies  of  sediment  transport 
processes.  For  this  reason,  it  is  important  that  wave  measurements  be  interpreted  as 
accurately  as  possible. 

Despite  the  need  for  good  data  on  the  kinematics  of  waves,  it  is  impractical  to 
measure  every  parameter  of  interest.  In  a  given  situation,  one  might  need  to  know  the 
velocities,  accelerations,  and  pressure  fluctuations  throughout  the  depth  at  a  given 
location.  While,  in  the  laboratory,  it  might  be  possible  to  measure  these  quantities 
at  many  locations,  it  would  be  very  expensive,  and  it  is  totally  impractical  in  the 
field.  A  pragmatic  approach  is  to  measure  just  a  few  quantities,  and  to  use  a  wave 
theory  to  compute  the  balance  of  the  parameters.  This  approach  can  give  a  complete 
description  of  the  kinematics  from  only  a  few  measurements.  Unfortunately,  the 
accuracy  of  this  approach  is  limited  by  the  accuracy  of  the  adopted  theory,  as  well  as 
the  accuracy  of  the  measurements  themselves. 
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A  number  of  wave  theories  have  been  adopted  in  an  effort  to  describe  the  kine¬ 
matics  of  waves.  The  most  accessible  and  frequently  used  of  these  is  Airy  wave  theory 
(also  known  as  linear  theory  or  first  order  Stokes  theory).  Airy  wave  theory  is  ea,sy 
to  use,  but  has  a  number  of  limitations.  The  source  of  these  limitations  is  the  sim¬ 
plifications  of  the  governing  equations  of  gravity  waves  that  are  made  to  linearize  the 
equations,  allowing  a  straightforward  solution.  These  simplifications  are  made  in  the 
free  surface  boundary  conditions  and  are  justified  by  an  assumption  of  small  wave 
amplitude.  Simplifying  the  free  surface  boundary  conditions  reduces  the  accuracy 
of  the  predicted  kinematics.  Unfortunately,  this  compromise  is  at  the  free  surface, 
which  is  the  location  of  the  greatest  fluid  velocities  and  accelerations  in  waves,  and 
thus  frequently  the  most  important  to  the  forcing  of  structures.  The  assumption  of 
small  amplitude  also  renders  the  theory  inadequate  for  large  waves,  exactly  those  of 
greatest  interest  to  coastal  engineers. 

In  order  to  address  these  limitations,  a  number  of  high  order  steady  wave  theo¬ 
ries  have  been  developed.  Commonly  in  use  are  Stokes,  Cnoidal,  and  Fourier  wave 
theories.  For  a  review  of  these,  see  Fenton  (1990).  In  general,  Stokes  methods  are 
successful  in  deep  water,  Cnoidal  methods  in  shallow  water,  and  Fourier  methods  in 
all  depths  of  water.  Within  their  limitations,  all  three  of  these  methods  provide  ex¬ 
cellent  predictions  of  the  kinematics  of  steady  waves,  but  are  not  directly  applicable 
to  the  irregular  waves  commonly  found  in  the  field.  With  the  possible  exception  of 
swell  conditions  on  a  very  mild-sloped  bottom,  waves  in  the  sea  arc  neither  steady, 
unidirectional  nor  monochromatic. 

Methods  for  the  interpretation  of  measurements  of  real  sea  states  rely  on  Airy  wave 
theory  even  more  heavily  than  do  steady  wave  methods.  With  modern  computer  sy.s- 
tems,  even  the  most  complicated  of  high  order  steady  wave  theories  is  quite  accessible. 
However,  such  higher  order  theories  are  not  directly  applicable  to  multi-directional  or 
multi-chromatic  waves.  The  linearity  of  Airy  theory  allows  superposition,  in  which 
any  combination  of  waves  of  different  frequencies  or  directions  can  be  combined  to 
form  a  solution.  Superposition  allows  real  sea  states  to  be  easily  characterized  by 
frequency  and  direction  spectra.  While  accessible,  this  method  results  in  solutions 
for  the  kinematics  of  the  waves  that  do  not  satisfy  the  full  free  surface  boundary 


"wave" 


Figure  1.1:  Segment  of  record  considered  for  global  approximations 


conditions,  and  thus  do  not  take  into  account  the  interactions  between  the  individual 
waves. 


1.2  Methods  for  the  Interpretation  of 
Irregular  Wave  Records 

Methods  used  for  the  interpretation  of  irregular  wave  records  fall  into  two  general 
categories:  global  and  local  approximations.  Global  methods  seek  a  solution  that 
matches  an  entire  measured  record,  or  a  single  complete  measured  wave,  from  trough 
to  following  trough,  or  zero  crossing  to  zero  crossing  (Fig.  1.1).  These  methods  apply 
the  same  frequency  and  wave  number  (or  set  of  frequencies  and  wave  numbers)  for 
all  z  (vertical  variation)  and  t  (time). 

Local  methods,  on  the  other  hand,  seek  an  approximation  to  each  small  local 
segment  of  a  measured  wave.  In  these  methods,  the  frequency  and  wave  number  still 
^PPly  for  all  z,  but  are  allowed  to  vary  with  time,  providing  a  separate  solution  in 
each  small  window  in  time  (Fig.  1.2). 
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Figure  1.2:  Segment  of  record  considered  for  local  approximations 


1.2.1  Global  Methods 

Spectral  Methods 

The  most  commonly  used  global  method  for  the  analysis  of  irregular  waves  is 
spectral  analysis,  coupled  with  superposition  of  linear  waves.  Spectral  methods  based 
on  a  single  point  measurement  seek  to  define  the  sea  state  with  a  variance  spectrum, 
E{u))  which  specifies  the  contribution  to  the  variance  of  the  water  surface  at  each 
frequency,  u).  In  practice,  the  spectrum  is  defined  in  discrete  form,  with  the  water 
surface  described  by  the  superposition  of  many  linear  waves, 

M 

ri{x,y,t)  =  ^  Cm  cos  {k^  (xcos^„  +  t/sin0„)  -  Umt  +  Om)  (1-1) 

m=l 

where  t]  is  the  elevation  of  the  water  surface,  and  km  are  the  frequency  and  wave 
number  of  the  mth  wave,  $m  is  the  direction  of  propagation  of  the  mth  w'ave  and 
is  the  phase  of  the  mth  wave  at  the  origin  of  the  coordinate  system. 

The  amplitudes  are  found  from  the  discrete  variance  spectrum, 

Qm  =  y/2E{uJm)^t^  (1.2) 

where  Au  is  the  frequency  spacing  of  the  discrete  spectrum.  The  frequency  and  wave 
number  of  each  of  the  M  waves  are  related  by  the  linear  dispersion  relation. 


(1.3) 
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where  g  is  the  acceleration  of  gravity,  and  h  is  the  mean  water  depth.  The  dispersion 
relation  defines  the  phase  speed  of  each  component,  allowing  each  separate  component 
to  move  independently,  without  interaction  with  one  another.  Any  Eulerian  current 
is  not  included  in  this  relation,  and  is  generally  ignored  in  spectral  methods. 

The  discrete  variance  spectrum,  E,  can  be  computed  from  a  water  surface  or  sub¬ 
surface  pressure  record  through  the  use  of  variations  of  the  Fast  Fourier  Transform 
(Dean  and  Dalrymple  1991;  Newland  1993).  A  point  measurement  provides  no  in¬ 
formation  about  the  direction  of  propagation  of  the  waves,  9^.  It  is  usually  assumed 
that  all  the  waves  are  propagating  in  the  same  direction.  Once  the  amplitude,  fre¬ 
quency,  wave  number  and  phase  of  each  individual  wave  are  defined,  the  kinematics 
and  dynamics  of  the  wave  field  can  be  computed  by  superimposing  the  kinematics  of 
each  individual  wave,  as  predicted  by  linear  wave  theory. 


Directional  Spectra  When  measurements  are  taken  by  an  array  of  instruments, 
the  directional  nature  of  the  sea  can  be  described  by  a  directional  spectrum,  S{uj,0). 
In  this  case,  the  water  surface  is  represented  by  a  large  number  of  linear  waves  of 
different  frequencies  and  directions: 


EE-.  n  COS  COS  0^  y  sin  0^^  —  QJm,n) 


where: 


^2S{um,en)Au;Ae 


and  Au;  and  A9  are  the  sample  spacing  of  the  discrete  spectrum  in  frequency  and 
direction  space.  u;„,  and  km  are  once  again  related  by  the  linear  dispersion  relation 
(Eq.  1.3).  The  directional  spectrum  is  usually  broken  down  into  two  parts,  the  one 
dimensional  variance  spectrum,  E{uj),  and  the  directional  spreading  function  (DSF), 
DM: 


S{uj,$)  =  E{u)D{u;,0) 


where  E{u>)  is  the  variance  spectrum  defined  above,  and 


D{uj,  $)d6du)  =  1 


(1.6) 

(1.7) 
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There  is  a  great  deal  of  literature  about  how  to  best  determine  the  DSF  from  a  variety 
of  arrangements  of  measurements.  Current  practice  is  reviewed  in  Mansard  (1997). 
Unfortunately,  all  of  these  methods  are  limited  in  their  ability  to  define  accurately 
the  DSF. 

The  time  series  measurements  used  to  determine  the  variance  spectrum  commonly 
include  in  excess  of  1000  points  in  time.  This  much  data  in  the  time  domain  allows 
a  very  high  resolution  computation  of  the  spectrum  of  a  stationary  process  in  the 
frequency  domain.  In  the  spatial  domain,  by  contrast,  there  are  only  as  many  data 
points  as  there  are  instruments  in  the  particular  measurement  array.  This  may  be 
as  few  as  three,  and  perhaps  as  many  as  a  dozen,  but  it  is  too  expensive  to  use 
many  more  than  that.  As  a  result,  there  is  little  information  to  define  the  DSF.  For 
example,  when  data  from  an  array  of  three  instruments  is  analyzed  by  the  standard 
linear  method,  the  DSF  can  be  specified  by  only  five  independent  coefficients  (Dean 
and  Dalrymple  1991).  While  these  few  coefficients  may  serve  well  for  determining 
integral  properties,  such  as  the  mean  direction  and  radiation  stress,  it  is  not  enough 
information  to  accurately  specify  the  complete  kinematics. 

Another  difficulty  arises  w'hen  determining  the  spectrum  from  measurements  other 
than  wave  staffs.  Most  often,  subsurface  pressure  gauges  or  a  combination  of  pressure 
gauges  and  orthogonal  velocity  gauges  are  used.  In  this  case,  the  measured  quantity 
must  be  related  by  a  transfer  function  to  the  equivalent  water  surface.  The  transfer 
function  is  most  commonly  determined  from  linear  wave  theory.  This  use  of  linear 
wave  theory  once  again  contributes  to  errors  in  situations  where  linear  theory  is  not 
entirely  appropriate.  Bishop  and  Donelan  (1987)  discuss  the  difficulties  in  determin¬ 
ing  the  transfer  function  from  a  subsurface  pressure  measurement  to  the  equivalent 
water  surface. 

The  commonly  used  methods  for  determining  the  DSF  are  statistical  methods, 
in  that  they  rely  on  the  assumption  that  the  phases  of  the  individual  components 
are  randomly  distributed.  This  assumption  allows  the  DSF  to  be  computed  by  dis¬ 
regarding  the  phaise  information.  Unfortunately,  without  the  phase  information,  it  is 
impossible  to  reconstruct  the  detailed  kinematics,  only  statistical  descriptions  can  be 
formulated. 
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Even  if  a  frequency  or  frequency-direction  spectrum  has  been  accurately  deter¬ 
mined,  these  methods  have  a  number  of  shortcomings  when  used  to  predict  the  kine¬ 
matics  of  irregular  waves.  Shortcomings  include  the  inaccuracies  inherent  in  linear 
wave  theory,  particularly  in  shallow  water  and  with  large  waves. 

An  additional  difficulty  arises  from  the  superposition  of  many  waves.  This  problem 
is  known  as  high  frequency  contamination  of  the  crest  kinematics  (Forristall  1985;  Lo 
and  Dean  1986;  Bishop  and  Donelan  1987).  Fundamentally,  the  difficulties  arise 
from  the  approximations  made  by  linear  wave  theory  in  the  free  surface  boundary 
conditions.  In  linear  theory  the  free  surface  boundary  conditions  are  applied  at  the 
mean  water  level,  and  thus  predictions  made  above  that  level  are  strictly  out  of 
the  solution  domain.  If  the  full  free  surface  boundary  conditions  are  not  satisfied, 
the  resulting  predictions  will  be  inaccurate,  particularly  near  the  free  surface.  In 
particular,  the  hyperbolic  function  quotients  that  define  the  vertical  variation  of  the 
kinematics  become  very  large  in  the  region  above  the  MWL  for  the  high  frequency 
(and  high  wave  number)  components.  This  results  in  substantial  high  frequency 
fluctuations  in  the  predicted  kinematics  near  the  crest. 

In  an  attempt  to  reduce  the  inaccuracies  in  linear  superposition’s  predictions  of  the 
near  surface  kinematics,  empirical  modifications  to  Airy  theory  have  been  adopted 
(Wheeler  1969;  Lo  and  Dean  1986).  This  method,  known  as  Wheeler  stretching, 
locally  adjusts  the  vertical  dimension  to  prevent  the  evaluation  of  the  hyperbolic 
quotients  from  being  evaluated  above  the  MWL.  The  result  is  a  hybrid  global-local 
method  in  that  the  frequency,  wave  numbers,  and  amplitude  of  each  wave  are  deter¬ 
mined  globally,  but  the  coordinate  system  is  defined  locally,  varying  with  time. 

The  stretching  method  produces  predictions  of  crest  kinematics  that  seem  to 
match  measured  data  better  than  simple  linear  superposition,  but  it  no  longer  satis¬ 
fies  either  the  Laplace  equation  (mass  conservation),  or  the  full  free  surface  boundary 
conditions. 

Another  attempt  to  improve  on  the  accuracy  of  determining  the  kinematics  from 
directional  spectra  is  the  recent  work  by  Prislin  et  al.  (1997),  and  Prislin  and  Zhang 
(1997).  This  method  seeks  to  reconstruct  wave  kinematics  from  a  directional  spec¬ 
trum  using  a  second  order  Stokes-type  interacting  wave  theory.  The  spectrum  is 
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decomposed  into  a  set  of  individual  free  waves,  with  from  one  to  five  directional  free 
waves  per  frequency.  The  effects  of  the  second  order  interactions  are  subtracted  from 
the  measured  record  to  determine  the  amplitudes  and  phases  of  the  individual  free 
waves  through  an  iterative  procedure.  This  results  in  an  expression  for  the  potential 
function  and  water  surface  that  includes  a  full  set  of  many  free  waves,  and  the  corre¬ 
sponding  second  order  bound  modes.  The  method  succeeds  in  globally  reproducing 
the  measured  kinematics  in  the  given  deep  water  field  records  quite  well,  although 
the  largest  errors  are  in  the  vicinity  of  the  crest,  where  velocities  ai(  greatest,  and 
accuracy  is  most  important. 

Another  limitation  of  the  Prislin  and  Zhang  method  is  that  the  nonlinear  interac¬ 
tions  are  computed  through  the  use  of  a  Stokes-type  perturbation  expansion  in  wave 
steepness.  Based  on  experience  with  steady  waves,  a  second  order  expansion  of  this 
type  is  likely  to  be  adequate  in  deep  water,  but  if  the  method  were  to  be  applied 
in  transitional,  and  especially  in  shallow  water,  a  much  higher  order  representation 
would  be  necessary  (Fenton  1990).  While  it  is  theoretically  possible  to  extend  the 
method  in  this  manner,  it  would  increase  the  magnitude  of  the  computation  substan¬ 
tially,  perhaps  prohibitively.  Representing  an  entire  record  in  the  global  sense  requires 
many  interacting  free  waves.  Considering  their  interactions  at  high  order  would  pro¬ 
duce  huge  numbers  of  interacting  terms  that  must  be  considered.  An  alternative, 
local,  approach  would  need  to  consider  far  fewer  interacting  waves. 

Other  global  methods  rely  on  zero  crossing  analysis  to  identify  particular  waves 
that  are  then  analyzed  b\'  using  steady  wave  theory  for  a  wave  of  the  same  height  and 
period.  This  approach  can  provide  an  order  of  magnitude  estimate  for  the  kinematics, 
but  does  not  take  into  account  the  detail  of  the  record,  and  thus  can  not  be  expected 
to  consistently  provide  better  than  order  of  magnitude  accuracy. 

To  include  the  detail  of  the  record  in  wave  by  wave  analysis,  Dean  (1965)  adapted 
his  numerical  stream  function  method  to  irregular  waves,  seeking  a  Fourier  expansion 
for  the  stream  function  that  includes  both  cos  j{kx  —  ivt)  and  smj{kT  —  ut)  terms. 
The  method  optimizes  the  Fourier  amplitudes  to  best  solve  the  free  surface  boundary 
conditions  at  the  measured  water  surface  elevations  from  trough  to  following  trough. 
While  this  approach  takes  into  account  the  detail  of  the  record,  the  method  includes 
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Figure  1.3:  Coordinate  system  for  Lambrakos-Baldock-Swan  method. 


only  a  single  free  mode,  with  all  other  included  components  being  bound  modes 
traveling  with  the  wave  at  the  same  phase  speed,  thus  representing  an  asymmetric 
wave  of  permanent  form.  A  solution  that  does  not  allow  a  change  in  form  is  unlikely 
to  accurately  capture  waves  in  deep  water,  where  frequency  dispersion  can  lead  to 
transient  extreme  waves,  or,  indeed,  in  shallow  water  where  shoaling  effects  cause  a 
change  in  form  as  the  waves  progress. 

Seeking  to  improve  on  global  methods,  Lambrakos  (1981)  developed  a  method  for 
determining  the  kinematics  of  two  dimensional  irregular  waves  that  includes  many 
free  modes,  and  thus  unsteady  motion.  Baldock  and  Swan  later  refined  Lambrakos’ 
method,  applying  it  specifically  to  large  transient  waves,  first  in  deep  water  (Baldock 
and  Swan  1994),  and  then  in  shallow  water  (Baldock  and  Swan  1996).  Baldock 
and  Swan’s  method  adopts  the  following  potential  function  that  is  a  double  Fourier 
expansion  in  space  and  time: 

N  M 

<^(x,  cosh  (nkz)  {An,m  COS  {nkx  -  mujt)  +  Bn,m  sin  {nkx  -  mut)) 

n=l  m=:l 

(1.8) 

where  k,  u),  global  constants,  and  ^  =  0  at  the  bed  (Fig.  1.3). 

In  this  potential  function  each  frequency  component  (mu;)  has  a  corresponding  set 
of  wavelengths  {nk),  allowing  each  component  to  travel  at  distinct  phase  speeds. 


10 


This  approach  can  accommodate  both  steady  and  unsteady  wave  forms  by  including 
both  a  set  of  free  waves  and  a  corresponding  set  of  bound  modes.  Periodicity  over 
time  and  space  is  assumed.  The  record  is  not  likely  to  be  periodic,  so  the  duration 
of  the  record  is  aissigned  as  the  fundamental  period,  (27r/uj),  and  the  fundamental 
wavelength,  {27^ jk)  is  found  as  part  of  the  solution.  The  method  uses  a  nonlinear 
optimization  to  find  the  coefficients,  and  Bn.m-,  the  fundamental  wavenumber, 
k,  and  the  Bernoulli  constant  that  produce  a  solution  that  globally  satisfies  the  full 
free  surface  boundary  conditions. 

In  order  to  accommodate  the  unsteadiness  of  the  wave  profile,  the  evolution  of 
the  wave  in  space  must  be  known.  This  is  accomplished  by  applying  the  free  surface 
boundary  conditions  at  many  locations  in  time  and  space.  Since  the  elevation  of  the 
water  surface  usually  is  measured  only  at  a  single  spatial  location,  it  is  predicted  at 
the  other  locations  as  part  of  the  solution.  The  potential  function  (Eq.  1.8)  exactly 
satisfies  the  bottom  boundary  condition  and  mass  conservation.  In  order  to  find  the 
solution  that  best  fits  the  free  surface  boundary  conditions,  the  method  seeks  a  set 
of  coefficients  that  minimizes  the  sum  of  squares  error  in  both  of  the  free  surface 
boundary  conditions  over  a  grid  of  nodes  in  time  and  space. 

Baldock  and  Swan  identified  a  difficulty  in  this  basic  method,  as  the  squared  error 
was  equally  considered  at  all  of  the  nodes,  most  of  which  were  at  locations  in  which 
the  water  surface  elevation  was  also  unknown.  This  resulted  in  solutions  in  which  the 
error  in  the  boundary  conditions  was  greatest  at  the  location  of  the  measurement. 
The  difficulty  was  mitigated  by  a  weighting  function,  multiplying  the  dynamic  free 
surface  boundary  condition  errors  at  the  measured  location  by  a  factor  of  50  in  the 
sum  of  squares  calculation,  to  force  the  solution  to  match  well  at  that  point. 

Comparisons  of  their  results  with  laboratory  data  were  quite  good,  but  the  method 
has  a  number  of  limitations.  These  include  a  huge  matrix  of  unknown  coefficients  that 
must  be  found  simultaneously  (2(A/ + 1)  unknowns,  with  typical  values  of  M  and  N 
of  18,  for  a  total  of  650  unknowns),  and  the  method  must  solve  for  the  water  surface  far 
from  the  measurement  location,  removing  the  focus  from  the  actual  measured  data. 
The  method  makes  no  assumptions  about  the  steadiness  of  the  wave  field;  however,  it 
extends  the  horizontal  bottom  assumption  and  introduces  unnecessary  complication 
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in  its  attempt  to  capture  an  entire  segment  of  an  unsteady  wave  field  with  a  single 
expression.  This  complication  would  become  at  least  an  order  of  magnitude  greater 
if  the  method  were  extended  to  include  the  second  horizontal  dimension. 


1.2.2  Local  Methods 

Nielsen  (1986,  1989)  introduced  a  method  that  uses  a  local  frequency  and  linear 
wave  theory  to  find  the  location  of  the  water  surface  from  a  pressure  record.  In  this 
method,  the  wave  in  the  region  of  each  measured  point  is  considered  to  be  a  small 
segment  of  a  linear  wave,  so  that  the  pressure  record  is  locally  represented  by  a  sine 
function. 


so  that: 


p  =  An  sin  (uJntn  “  ^n) 


(1.9) 


OJn 


1  d'^p 

p  df^ 


(1.10) 


Estimating  the  curvature  by  finite  differences,  the  local  frequency  becomes: 


iOn  = 


~Pn—l  ”t"  2pn  Pn+1 


Pn^t^ 


(1.11) 


where  At  is  the  time  spacing  between  points.  Alternatively,  u  can  be  computed  from 
trigonometric  identities. 


LOn  =  -  cos' 


-1  /  Pn-l  +  Pn+1  \ 

V  2p„  ) 


(1.12) 


which  is  exact  for  a  sinusoid. 

Once  the  frequency  has  been  determined,  the  wave  number  is  computed  from  the 
linear  dispersion  relation,  Eq.  1.3.  The  water  surface  elevation  can  then  be  computed 
with  the  linear  pressure  response  function: 


Pn  COsh{knh) 

nn  = - TT^ - T  1.13) 

pg  cos\i{knZp)  ''  ’ 

where  Zp  is  the  elevation  of  the  pressure  measurement  above  the  bed,  and  h  is  the 
mean  water  depth  (As  the  method  is  a  local  method  for  the  interpretation  of  a  point 
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measurement,  the  local  mean  water  depth  is  used,  rather  than  the  global  still  water 
depth).  This  may  be  adapted  somewhat  by  using  a  variation  of  Wheeler  stretching: 

^  Pr.  ^  +  S)) 

pg  cosh(A-ni:p) 

Nielsen  also  presented  an  alternative  approach,  which  uses  a  semi-empirical  trans¬ 
fer  function  derived  from  Fourier  wave  theory  to  compute  the  water  surface  elevation 
from  the  measured  pressure  and  local  curv'ature  of  the  pressure  record.  While  quite 
accurate  at  reproducing  the  water  surface,  neither  method  supplies  the  kinematics, 
nor  do  they  satisfy  either  mass  conservation  and  or  the  free  surface  boundary  con¬ 
ditions.  Despite  these  limitations,  the  efficacy  of  these  methods  demonstrates  the 
potential  for  the  local  approach. 

To  interpret  bottom  pressure  measurements  in  the  context  of  the  kinematics  while 
preserving  the  full  governing  equations,  Fenton  (1986)  presented  a  method  that  em¬ 
ploys  a  local  polynomial  approximation  to  the  complex  potential  function.  In  Fenton’s 
method,  the  potential  function  and  water  surface  are  represented  by  separate  poly¬ 
nomials  in  each  small  window  in  time: 

M 

<f>{T  -  ct,  y)  -f  -  ct,  y)  =  ~  (1-15) 

M 

j=0 

where  2  =  x  -t-  ty,  y  =  0  at  the  bed,  t)  is  the  water  surface,  c  is  the  wave  celerity, 
and  the  aj  and  bj  are  real.  The  wave  is  assumed  to  propagate  at  speed,  c,  without 
change  in  form.  While  steadiness  is  not  a  valid  assumption  in  the  global  sense,  it  is 
only  applied  locally,  within  a  small  window  in  time. 

The  method  solves  for  the  coefficients,  aj  and  6_,,  and  the  wave  celerity,  c,  that 
satisfy  the  full  nonlinear  free  surface  boundary  conditions  and  fit  the  measured  pres¬ 
sure  record.  This  approach  provides  the  complete  kinematics  and  satisfies  the  full 
governing  equations.  Based  on  a  polynomial  variation  with  depth,  it  works  well  in 
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shallow  water,  where  Cnoidal  methods  with  polynomial  vertical  variation  are  theo- 
reticall}'  appropriate.  At  the  same  time,  it  may  not  be  applicable  in  transitional  or 
deep  water,  where  Stokes  type  methods,  based  on  an  exponential  variation  in  the 
vertical,  are  theoretically  more  appropriate.  In  a  later  paper,  Fenton  and  Christian 
(1989)  presented  a  simplified  version  of  the  method  that  required  less  algebraic  ma¬ 
nipulation,  and  was  still  found  to  be  effective  in  shallow  water.  In  this  case,  the  local 
wave  celerity  was  assumed  to  be  that  given  by  long  wave  theory: 

c=v^  (1.17) 

While  this  is  a  reasonable  approximation  for  long  waves  in  shallow  water,  it  was 
not  appropriate  for  shorter  waves  in  deeper  water,  as  might  be  expected  from  the 
polynomial  form  and  the  long  wave  celerity. 

To  find  a  method  that  could  work  in  any  depth  of  water,  Sobey  (1992)  developed 
the  Local  Fourier  Method  for  Irregular  waves  (LFI).  This  approach  employs  a  po¬ 
tential  function  represented  by  a  low  order  Fourier  expansion  in  a  small  window  in 
time.  It  is  a  method  derived  for  the  analysis  of  a  point  water  surface  trace.  Local 
frequency,  wave  number,  and  the  Fourier  coefficients  are  sought  that  fit  the  measured 
record  and  the  full  free  surface  boundary  conditions.  The  LFI  method  provides  the 
complete  kinematics,  satisfies  the  full  governing  equations,  and  is  successful  in  all 
depths  of  water.  A  more  complete  description  of  this  method  follows  in  Chapter  2. 
Sobey’s  method  shows  a  great  deal  of  promise,  but  was  only  applied  to  the  analysis 
of  a  water  surface  measurement  at  a  single  location. 

Pressure  sensors  are  easy  to  deploy,  and  thus  are  often  used  to  measure  waves, 
particularly  in  shallow  and  transitional  depth  water.  The  LFI  method  is  extended 
in  this  dissertation  to  the  interpretation  of  wave  measurements  from  a  point  pressure 
sensor  in  Chapter  3. 

Point  measurements  provide  no  information  about  the  directionality  of  the  mea¬ 
sured  wave  field.  Wind  seas  are  not  uni-directional,  and  knowledge  of  directionality 
has  been  shown  to  be  very  important  in  the  prediction  of  the  kinematics  and  forcing 
of  structures  by  real  seas  (Dean  1977;  Forristall  et  al.  1978). 

Arrays  of  measurements  are  frequently  used  in  order  to  capture  the  directional 
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nature  of  the  wave  field.  The  LFI  method  has  been  further  extended  in  this  disserta¬ 
tion  to  the  interpretation  of  arrays  of  water  surface  measurements  in  Chapter  5,  and 
arrays  of  pressure  sensors  in  Chapter  6. 
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Chapter  2 

Two  Dimensional  LFI  Theory 


It  is  common  to  deploy  a  single  measurement  device  to  gather  information  about 
the  wave  field  at  a  given  location.  Measurement  at  a  single  point  in  space  does  not 
give  any  information  about  the  directionality  of  the  sea  state,  but  if  the  wave  field 
is  assumed  to  be  locally  two  dimensional,  a  reasonable  description  of  the  kinematics 
can  be  established.  This  chapter  outlines  the  LFI  method  as  it  is  applied  to  the 
interpretation  of  a  time  series  from  measurements  taken  at  a  single  location  in  space, 
such  as  recorded  by  a  pressure  gauge  or  wave  staff. 


2.1  Problem  Formulation 


The  problem  formulation  for  two  dimensional  irregular  waves  has  much  in  common 
with  classical  steady  wave  theory.  The  flow  is  assumed  to  be  incompressible  and 
irrotational.  The  kinematics  can  therefore  be  represented  by  a  potential  function, 
(f>{x,z,t),  where 


and  u  and  w  are  the  horizontal  and  vertical  velocities,  respectively. 


(2.1) 
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Figure  2.1:  Coordinate  system  used  for  the  two  dimensional  method 


Field  Equation  The  field  equation  is  mass  conserv^ation  for  irrotational  flow,  in 
the  form  of  the  Laplace  equation: 


vV 


0J-2 


(2.2) 


Boundary  Conditions  The  boundary  conditions  are  the  bottom  boundary  condi¬ 
tion  for  a  locally  horizontal  bed  (BBC), 

u>  =  ^  =  0  at  z  =  —h  (2.3) 

the  kinematic  free  surface  boundary  condition  (KFSBC), 

d?!  dr] 

and  the  dynamic  free  surface  boundary  condition  (DFSBC), 

^  +  +  gr]-B  =  0  at  z  =  r}  (2.5) 

where  t?  is  the  elevation  of  the  free  surface,  g  is  the  acceleration  of  gravity,  and  B  is 
the  Bernoulli  constant. 

The  kinematic  free  surface  boundary  condition  (Eq.  2.4)  includes  both  the  time 
and  spatial  gradients  of  the  water  surface.  When  working  with  a  subsurface  gauge. 
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the  location  of  the  water  surface  is  unknown  and  part  of  the  solution.  When  working 
with  data  from  a  water  surface  probe,  the  location  of  the  water  surface  at  discrete 
points  in  time  is  known,  but  the  gradients  in  time  and  space  are  not.  In  order 
to  accommodate  this  lack  of  data,  the  kinematic  free  surface  boundary  condition 
is  transformed  to  eliminate  the  gradients  of  the  water  surface.  Following  Longuet- 
Higgins  (1962),  the  kinematic  condition  is  subtracted  from  the  the  dynamic  condition 
differentiated  following  the  motion: 


MKFSBC  =  -5r(KFSBC)  +  -^(DFSBC) 

JL/  L 

Resulting  in  a  modified  kinematic  free  surface  boundary  condition: 

du  dw. 

^  +  gw  +  2(u-^+v^-) 

2du  dw 

+  u  —  +  uw— 
ox  ox 


,  du 

+  uw— - 1-  w 

Oz 


,dw 


at 


z  =  rj 


(2.6) 


This  new  form  of  the  boundary  condition  does  not  include  the  gradients  of  the  water 
surface.  Applying  both  the  modified  kinematic  free  surface  boundary  condition  and 
the  dynamic  free  surface  boundary  condition  completes  the  formulation. 


Observational  Equations  In  steady  wave  theory,  periodic  boundary  conditions 
are  also  imposed,  forcing  the  solution  to  be  periodic  in  both  space  and  time.  For 
irregular  waves,  however,  the  periodicity  is  not  known.  Rather,  the  local  solution 
is  defined  by  a  local  segment  of  a  measured  record  within  a  small  window  in  time, 
together  with  the  field  equation  and  the  bottom  and  free  surface  boundary  conditions. 

In  order  to  define  a  solution  that  fits  the  measured  record,  observational  equations 
are  identified.  These  equations  will  be  different  depending  on  which  quantity  has 
been  measured.  In  the  case  of  a  water  surface  measurement,  they  are  the  free  surface 
boundary  conditions,  applied  at  the  measured  elevations  at  a  number  of  points  in 
time  throughout  the  window  (Sobey  1992).  In  the  case  of  a  pressure  measurement, 
they  are  the  Bernoulli  equation,  applied  at  the  elevation  of  the  measurement,  and 
also  at  a  number  of  points  in  time  within  the  window  (see  Chapter  3). 
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Solution  in  Each  Local  Window  A  form  for  the  potential  function  in  each  local 
window  is  based  on  Fourier  steady  wave  theory.  This  is  the  potential  function  used 
by  Sobey  (1992): 


<i>[x,z,t)  =  Ux  +  ^  Aj 
;=i 


cosh  jk{h  + 
cosh  jkh 


s\r\j{kx  —  u:t) 


(2.7) 


where  U  and  h  are  the  known  depth-uniform  Eulerian  current  and  mean  w’ater  depth 
(see  section  3.3.1  for  a  discussion  of  these  important  parameters),  J  is  the  truncation 
order  of  the  Fourier  series,  Aj  are  the  Fourier  coefficients,  and  w  and  k  are  the 
local  fundamental  frequency  and  wave  number.  The  above  potential  function  exactly 
satisfies  mass  conservation  and  the  BBC.  This  form  for  the  potential  function  is 
periodic  in  space  and  time,  however  the  periodicities  are  not  defined  apriori,  but 
found  to  fit  the  record,  defining  a  local  frequency  and  wave  number. 

The  measured  record  is  broken  down  into  individual  segments,  each  in  a  separate 
window  in  time.  In  each  window,  a  different  set  of  the  parameters  w,  fc,  kx,  and  Aj 
are  found  to  fit  the  segment  of  the  measured  record.  This  represents  that  segment  as 
a  piece  of  a  larger,  periodic  w'ave.  The  entire  record  is  then  represented  by  separate 
potential  functions,  each  applied  to  a  particular  window  in  time. 


2.1.1  Dynamics 

The  potential  function  provides  the  complete  kinematics,  and  the  dynamics  are 
found  through  the  unsteady  Bernoulli  equation: 

+  -(u^  +  u'^)  -h  gz  -  B  =  Q  (2.8) 

where  p  is  the  mass  density  of  the  water,  and  p  is  total  pressure.  Total  pressure  below 
the  water  surface  can  be  broken  down  into  three  components.  Atmospheric  pressure 
(Pa)  is  the  pressure  of  the  atmosphere  at  the  w’ater  surface,  hydrostatic  pressure  (p/,) 
is  the  pressure  in  the  w'ater  column  due  to  gravity,  in  the  absence  of  motion.  Dynamic 
pressure  (pd)  is  the  component  due  to  the  motion  of  the  fluid.  Separating  these  three 
components  focuses  attention  on  the  wave  motion. 

p  =  Pa  +  P<f  +  P/i  where  ph  =  -pgz  for  z  <t]  (2.9) 
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When  atmospheric  pressure  is  defined  as  zero,  and  dynamic  pressure  is  substituted 
into  Eq.  2.8,  the  result  is: 

+  +  +  ^ -5  =  0  (2.10) 

The  Bernoulli  Constant 

In  a  potential  flow,  the  Bernoulli  constant  is  the  same  throughout  space  and 
time.  For  the  special  case  of  wave  motion,  the  value  of  the  Bernoulli  constant  can 
be  computed  if  the  kinematics  are  known  (Longuet-Higgins  1975).  The  Bernoulli 
equation  (2.8)  and  the  bottom  boundary  condition  are  applied  at  the  bottom: 

d<j>  1  2  Ph  — 

■5J-+2%+9^6  +  --B  =  0  (2.U) 

where  u;,,  pj,  and  are  the  velocity,  pressure,  and  elevation  at  the  sea  bed.  If  the 
flow  is  periodic,  a  time  average  over  a  period  results  in: 

l~o  Vk  — 

-ul  -\-  gzb  +  ^  -  B  =  Q  (2.12) 

where  the  over-bars  indicate  time  averaging. 

In  the  case  of  steady  periodic  wave  motion,  the  total  vertical  momentum  at  any 
horizontal  location  must  be  the  same  at  the  beginning  and  end  of  a  period.  In  order 
for  this  to  be  the  case,  the  vertical  momentum  averaged  over  a  period  must  be  a 
constant.  In  order  for  the  momentum  to  remain  constant,  the  time  averaged  net 
vertical  force  on  the  water  column  at  that  point  must  be  zero.  If  the  pressure  at  the 
water  surface  {pa)  is  taken  to  be  zero,  then  the  force  of  the  pressure  at  the  bed  must 
be  equal  to  the  gravitational  force  on  the  column  of  water,  so  that  the  mean  pressure 
on  the  bed  is  hydrostatic, 

Pb  =  P9{v-Zb)  (2.13) 

resulting  in  a  simple  and  exact  expression  for  the  Bernoulli  constant. 

B  =  Ji?  +  iuj 


(2.14) 
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With  2  defined  to  be  zero  at  the  mean  water  level  (rj)  and  Eq.  2.7  as  the  potential 
function,  B  becomes  (Sobey  1992): 


B  = 


y(  jkA,  y 
Y  \coshj kh  J 


(2.15) 


and  all  the  terms  in  Eq.  2.10  are  defined  by  the  potential  function,  allowing  the 
dynamic  pressure  to  be  computed  from  the  kinematics. 


2.2  Finding  the  Solution 

The  bulk  of  the  LEI  method  is  the  process  of  determining  appropriate  values  for 
each  of  the  parameters,  u>,  k,  kx,  and  Aj,  in  the  potential  function  (Eq.  2.7)  for  each 
window  in  time.  The  result  is  a  set  of  potential  functions  that  completely  describe  the 
kinematics  of  the  wave  field  in  the  region  of  the  measurement,  each  separate  window 
in  time  being  described  by  a  unique  and  independent  potential  function. 

In  the  case  of  subsurface  measurements,  the  location  of  the  water  surface  is  also 
required.  Determination  of  the  parameters  of  the  potential  function,  as  well  as  the 
location  of  the  water  surface,  is  accomplished  through  nonlinear  optimization  that 
matches  the  measured  record  and  minimizes  the  error  in  both  free  surface  boundary 
conditions.  Sobey  (1992)  applied  the  LEI  method  to  the  interpretation  of  a  single 
water  surface  trace.  The  following  chapter  applies  the  method  to  a  subsurface  pressure 
record. 
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Chapter  3 

LFI  Method  for  a  Point  Pressure 
Record 


Pressure  gauges  are  relatively  inexpensive  and  easy  to  deploy,  and  as  a  result  are 
frequently  used  in  the  field  to  measure  waves,  particularly  in  shallow  water.  The 
simplest  instrument  consists  of  a  single  subsurface  pressure  gauge.  While  a  single 
point  measurement  provides  no  directional  information,  a  reasonable  estimate  for  the 
kinematics  of  the  measured  waves  can  be  obtained.  This  chapter  describes  a  technique 
for  applying  the  two  dimensional  LFI  method  presented  in  the  previous  chapter  to 
the  analysis  of  a  point  subsurface  pressure  trace. 


3.1  Formulation  of  Solution 

The  flow  is  taken  to  be  two  dimensional,  with  the  governing  equations  described 
in  Chapter  2.  These  include  the  potential  function. 


M)  =  (3.,) 

the  modified  kinematic  free  surface  boundary  condition  (/^ ), 
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/"■  = 


d^d  du  dw 

2  du  dir 
+  w  —  +  ttw-^ 
dx  dx 


du 

4-  uw— — h  w 
dz 


,5u’ 


=  0 


at  z  =  -q 


(3.2) 


the  dynamic  free  surface  boundary  condition  (/^), 


/^  =  ?^  +  -\-gq-B  =  {i  at  z-q  (3.3) 

and  the  unsteady  Bernoulli  equation  (/®). 

/*  =  ^  +  5("’  +  >-'‘)  +  ^-B  =  0  (3.4) 

with  the  Bernoulli  constant  defined  as: 


The  unknown  parameter,  j,  appears  in  the  potential  function  onlj'  when  coupled 
with  the  parameter,  k.  It  is  convenient  to  solve  for  the  non-dimensional  parameter, 
Arx,  essentially  a  phase  parameter  in  the  potential  function. 

The  action  of  waves  is  greatest  near  the  surface,  so  that  the  fluid  velocities  and 
accelerations  are  largest  near  the  surface  and  decay  rapidly  with  depth,  particularly 
in  deep  water.  In  addition,  for  pressure  sensors  placed  at  or  near  the  bottom,  the 
vertical  component  of  the  velocity  is  damped  out  completely,  as  expressed  by  the 
bottom  boundary  condition.  The  problem,  then,  is  to  extrapolate  the  motion  of  the 
fluid  up  to  the  surface,  where  the  motion  is  greatest,  using  only  data  extracted  from 
the  subsurface,  where  the  motion  is  much  less. 
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Mathematically,  this  is  an  ill  posed  problem,  as  the  flow  is  governed  by  an  elliptic 
equation  (Laplace,  Eq.  2.2),  in  which  the  solution  is  determined  by  the  boundaries, 
but  there  are  only  data  at  a  single  location.  Some  of  these  difficulties  can  be  sur¬ 
mounted  through  the  application  of  the  free  surface  boundary  conditions.  While  there 
are  no  data  measured  near  the  surface,  the  free  surface  boundary  conditions  remain 
appropriate,  and  necessary  to  define  the  solution  at  that  boundary. 

When  the  LFI  method  was  applied  to  a  water  surface  trace  (Sobey  1992),  the 
location  of  the  the  water  surface  was  known,  and  the  boundary  conditions  could 
be  directly  applied  at  that  location.  Unfortunately,  when  working  from  a  subsurface 
pressure  record,  the  location  of  the  water  surface  is  not  known.  W^hile  the  free  surface 
boundary  conditions  are  well  defined,  the  fact  that  the  location  at  which  they  must 
be  applied  is  not  known  makes  the  problem  more  difficult.  This  is  the  complication 
that  leads  to  the  difficulties  in  finding  full  nonlinear  solutions  to  all  free  surface  flow 
problems.  In  order  to  apply  the  free  surface  boundary  conditions  when  working  with 
a  subsurface  record,  the  location  of  the  water  surface  must  be  found,  together  with 
the  potential  function,  in  each  window. 

In  order  to  locate  the  free  surface,  the  water  surface  is  defined  at  N  surface 
nodes  equally  spaced  in  time  throughout  the  window.  The  elevation  of  these  nodes 
is  unknown,  and  will  be  sought  as  part  of  the  solution.  Including  the  water  surface 
nodes  as  part  of  the  sought  solution  introduces  N  additional  unknown  parameters  for 
a  total  of  3  -f  J  -f-  unknowns  in  each  window  (k,  kx,  u,  Ai . . .  Aj,  rji . . .  r]N).  The 
free  surface  boundary  conditions,  Eq.  3.3  and  3.2,  are  applied  at  each  surface  node. 


3.2  Formulation  of  the  Optimization 

Finding  the  unknown  parameters  in  a  nonlinear  system  of  algebraic  equations  is 
known  as  nonlinear  optimization.  The  system  in  this  case  consists  of  the  nonlinear 
Bernoulli  equation  and  the  nonlinear  free  surface  boundary  conditions.  The  free 
surface  boundary  conditions  are  nonlinear  in  two  respects,  involving  second  order 
terms  in  dynamic  free  surface  boundary  condition  (the  terms  and  B  in  Eq.  3.3), 
and  second  and  third  order  terms  in  the  modified  kinematic  boundary  condition  (the 
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MWL 


Elevations  of  water  surface 
nodes  are  sought 


Bernoulli  equation  applied 
at  known  elevation,  zp 


Figure  3.1:  Schematic  of  system  of  equations  in  a  window 

and  terms  in  Eq.  3.2).  There  is  also  a  significant  nonlinearity  introduced  by 
the  application  of  the  boundary  conditions  at  the  unknown  and  varying  free  surface. 

Observational  Equations  The  given  form  for  the  potential  function  could  repre¬ 
sent  any  periodic  flow,  subject  to  the  bottom  boundary  condition.  The  FSBCs  define 
the  flow  as  a  gravity-constrained,  free  surface  flow.  The  observational  equations  are 
the  equations  in  the  system  that  force  the  solution  to  fit  the  given  record.  For  a 
subsurface  pressure  record,  this  is  the  Bernoulli  equation,  applied  at  the  location  of 
the  pressure  measurement.  The  required  number  of  independent  equations  are  es¬ 
tablished  by  applying  the  Bernoulli  equation  at  a  number  of  times  throughout  the 
window  considered.  The  error  in  the  Bernoulli  equation  is  the  difference  between  the 
measured  dynamic  pressure  and  that  computed  from  the  kinematics  defined  by  the 
potential  function.  The  solution  is  the  set  of  parameters  in  the  potential  function 
and  the  set  of  water  surface  nodes  that  produces  a  predicted  dynamic  pressure  that 
matches  the  measured  record,  while  simultaneously  satisfying  the  FSBCs. 

A  system  of  equations  is  specified  if  there  are  as  many  independent  equations  as 
unknown  parameters  in  the  system.  If  there  are  more  equations  than  unknowns,  the 
solution  can  be  defined  as  that  which  results  in  the  smallest  squared  errors  in  the 
equations.  This  least  squares  formulation  is  also  appropriate  for  a  uniquely  defined 
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system,  with  an  expected  error  of  zero  in  all  the  equations.  In  order  to  specify  the 
solution  to  the  LFI  formulation,  the  boundary  conditions  (f^  and  are  applied  at 
each  of  the  water  surface  nodes,  yielding  2N  equations,  and  the  Bernoulli  equation 
ifF)  is  applied  at  /  nodes  on  the  pressure  record  within  the  local  window  (Fig.  3.1). 
There  are  a  total  of  3  +  J  +  A'’  unknowns  {k,  kx,  u,  Ai ...  Aj,  and  r)i...  r]N).  The 
system  is  uniquely  specified  for  /  +  =  3  +  J  and  overspecified  for  /  +  A"  >  3  +  J, 

resulting  in  the  following  least  squares  optimization: 

N 

minimizeC>(X)  =^/f  + /f  (X;7y„,t„)2 

n=l 

/ 

+  (3.6) 

t=l 

X  =  {k,  kx,  uj,Ai...Aj,T]i...  ijn) 

where  are  the  locations  in  time  where  the  water  surface  nodes,  rjn,  are  sought. 
ti  are  the  times  within  the  window  where  the  Bernoulli  equation  is  applied,  Pd, 
are  the  measured  dynamic  pressures  at  ti,  and  Zp  is  the  elevation  of  the  pressure 
gauge.  Overspecification,  particularly  with  additional  nodes  on  the  pressure  record, 
is  advantageous  for  an  actual  record  as  a  way  to  minimize  the  effect  of  any  noise  in 
the  measurements.  Additional  nodes  on  the  pressure  record  also  serve  to  emphasize 
the  measured  data,  which  directly  define  the  local  kinematics.  The  details  of  the 
solution  to  this  system  of  equations  is  given  in  the  following  section. 


3.3  Computation  Methods 

The  LFI  method  can  be  broken  down  into  the  following  sequence  of  steps: 

1.  Pre-processing  of  record. 

(a)  Determine  estimate  for  level  of  noise  in  the  record. 

(b)  Determine  MWL  and  subtract  hydrostatic  pressure  from  record. 

(c)  Determine  estimate  for  magnitude  and  direction  of  Eulerian  current. 
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(d)  Define  a  continuous  record  from  the  discrete  observations  by  cubic  spline 
interpolation. 

(e)  Specify  spacing  of  output  locations. 

(f)  Compute  mean  zero  crossing  frequency. 

(g)  Non-dimensionalize  record  and  all  parameters. 

2.  Primary  values  of  numerical  solution  parameters  are  chosen. 

(a)  Window  width  (tq) 

(b)  Order  of  solution  {J) 

(c)  Number  of  sample  locations  on  record  within  each  window  (/) 

(d)  Number  of  water  surface  nodes  within  each  window  {N) 

3.  For  each  selected  output  location,  a  window  in  the  record  is  defined,  and  an 
LFI  solution  is  computed. 

(a)  Initial  guess  for  the  optimization  is  computed. 

(b)  Full  nonlinear  optimization  for  all  unknown  components  of  the  potential 
function  is  computed. 

(c)  Results  are  checked  for  spurious  solution. 

i.  If  no  solution,  or  a  spurious  solution,  is  found  the  solution  parameters 
are  adjusted,  and  the  optimization  repeated. 

ii.  If  a  good  solution  is  found,  progress  to  the  next  window. 

3.3.1  Pre-Processing  of  Record 

Accommodating  Measurement  Error 

Measurement  error  can  be  a  major  source  of  difficulty,  as  the  method  relies  on 
the  details  of  a  small  segment  of  the  record.  In  Sobey’s  work  (1992),  the  primitive 
kinematic  free  surface  boundary  condition  (Eq.  2.4)  was  applied  at  the  measured  water 
surface,  requiring  an  estimate  for  the  time  gradient  of  the  surface.  The  estimation 
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of  the  gradients  from  a  measured  record  is  very  sensitive  to  error  in  the  record,  and 
thus  Sobej^  found  it  necessary  to  smooth  field  and  laboratory  measurements  with  a 
simple  moving  average  filter. 

There  are  no  pressure  gradient  terms  in  the  Bernoulli  equation,  so  the  application 
of  the  LFI  method  to  a  measured  pressure  trace  should  be  less  sensitive  to  noise. 
Nevertheless,  measurement  error  remains  a  problem  in  applying  the  LFI  method  to  a 
pressure  record.  When  there  is  obvious  noise  in  the  record,  an  alternative  to  smooth¬ 
ing  is  to  substantially  overspecify  the  system  of  equations.  By  taking  many  samples 
on  the  pressure  record  in  each  window,  the  least  squares  optimization  may  accommo¬ 
date  the  errors  in  the  measurements.  This  is  preferable  to  smoothing,  as  no  smoothing 
algorithm  can  reproduce  lost  data,  but  will  rather  impose  some  apriori  assumption 
about  the  nature  of  the  record  on  the  results.  Accommodating  the  measurement  error 
with  many  samples  on  the  pressure  record  allows  the  least  squares  optimization  to 
find  the  best  fit  without  imposing  any  further  assumptions  on  the  results.  This  was 
the  approach  taken  with  the  presented  laboratory  data  (Section  3.5).  Records  gener¬ 
ated  analytically  by  Fourier  wave  theory  contained  no  error,  and  therefore  required 
no  special  accommodation. 

The  Mean  Water  Level 

The  mean  water  depth,  h,  is  included  in  the  potential  function  (Eq.  3.1).  This 
IS  an  unknown  value  when  the  only  data  available  is  a  pressure  record.  While  an 
approximate  depth  of  water  is  likely  to  be  known  at  a  given  deployment  site,  the 
water  level  fluctuates  due  to  astronomical  and  storm  tides.  These  processes  happen 
over  a  much  larger  time  scale  than  the  periods  of  wind  generated  waves.  An  estimate 
for  the  local  mean  water  level  (MWL)  can  be  determined  by  averaging  the  measured 
pressure  (after  subtracting  out  the  atmospheric  pressure)  over  a  period  of  time  that 
must  be  larger  than  a  typical  wave  period,  but  smaller  than  the  period  of  tidal 
fluctuations,  to  accommodate  changes  in  the  mean  water  level.  This  mean  pressure 
IS  the  local  hydrostatic  pressure,  which  is  used  to  compute  the  mean  water  level,  and 
is  subtracted  from  the  pressure  record  to  define  the  dynamic  pressure.  While  this  is 
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not  an  exact  procedure,  the  MWL  is  not  a  critical  parameter,  as  h  is  only  used  to 
locate  the  origin  of  coordinates  in  the  potential  function.  Note  that  the  result  is  a 
local  MWL,  which  could  be  quite  different  than  the  global  still  water  level  (SWL). 
The  mean  w'ater  level  is  the  more  appropriate  reference  frame,  as  the  method  is  local, 
and  only  considers  a  small  window  in  time  at  once.  The  elevations  of  the  water 
surface  nodes  are  found  independently,  without  any  requirement  that  the  resulting 
mean  water  level  be  at  the  origin.  Thus  it  is  not  necessary  that  the  origin  is  at  the 
exact  MWL,  only  that  the  it  is  near  the  surface.  In  shallow  water,  the  bed  could  be 
used  as  the  origin  (Fenton  1986),  but  this  would  not  work  well  in  deeper  water,  as 
the  area  of  primary  interest,  the  water  surface,  would  be  far  from  the  origin. 

The  Eulerian  Current 

The  depth  uniform  Eulerian  current,  t/,  appears  as  a  known  parameter  in  the 
potential  function,  Eq.  3.1.  Unfortunately,  a  measured  pressure  record  gives  no  in¬ 
dication  of  the  local  current.  Unlike  the  mean  water  level,  the  Eulerian  current  is  a 
critical  parameter  that  defines  the  propagation  medium,  and  varying  the  value  used 
for  the  current  will  have  a  considerable  effect  on  the  results  of  the  computation.  In 
fact,  as  pointed  out  by  Fenton  (1985a),  any  attempt  to  apply  a  wave  theory  without 
knowledge  of  the  local  current  will  be  inaccurate,  even  at  only  first  order. 

The  situation  is  not  hopeless,  however.  The  presence  of  the  Eulerian  current 
term  in  the  potential  function  draws  attention  to  the  fact  that  it  is  an  important 
parameter,  even  if  it  is  taken  to  be  zero.  This  should  encourage  investigators  planning 
field  experiments  to  establish  a  method  of  estimating  the  local  current.  This  may  be 
as  straightforward  as  placing  a  current  meter  nearby,  or  as  simple  as  using  tidal 
data  to  estimate  the  local  current.  Caution  must  be  taken  with  the  later  method, 
as  the  accuracy  of  the  estimate  may  not  be  consistent  with  the  use  of  high  order 
interpretation  methods. 

It  is  also  assumed  in  the  LFI  method  that  local  current  is  depth  uniform.  Although 
this  is  unlikely  to  be  strictly  true,  Kishida  and  Sobey  (1988)  demonstrated  that  for 
steady  waves,  a  current  profile  w’ith  a  realistic  level  of  vorticity  was  likely  to  yield  very 
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similar  results  as  a  depth  uniform  current  profile  for  the  wave  generated  kinematics, 
even  at  high  order.  It  is  expected  that  these  results  are  applicable  to  irregular  waves 
as  well,  and  thus  the  assumption  of  depth  uniform  current  is  reasonably  appropriate. 

Another  complication  to  be  considered  is  the  direction  of  any  local  current  with 
respect  to  the  propagation  direction  of  the  waves.  The  U  in  Eq.  3.1  is  the  component 
of  the  local  current  in  the  direction  of  wave  propagation.  In  the  case  of  a  single 
point  measurement,  no  information  is  available  to  indicate  the  wave  direction.  Even 
if  the  current  magnitude  and  direction  are  known  exactly,  without  an  estimate  for  the 
propagation  direction  of  the  waves,  an  accurate  estimate  for  U  is  not  available.  Once 
again,  further  information  could  be  used  to  help  resolve  this  difficulty.  Knowledge 
of  the  local  bathymetry  and  the  wave  climate  could  provide  this  information.  For 
example,  waves  near  a  beach  tend  to  propagate  almost  perpendicular  to  the  shore,  or 
if  there  are  directional  measurements  taken  nearby  available,  they  could  be  combined 
with  refraction  computations  to  provide  an  estimate  of  the  local  propagation  direction. 
There  must  be  some  information  available  about  both  the  local  Eulerian  current 
and  the  local  wave  propagation  direction  in  order  to  interpret  accurately  a  point 
measurement. 

The  pressure  records  used  for  the  examples  in  this  chapter  were  generated  either 
analytically  or  in  a  laboratory  flume.  In  both  cases,  the  Eulerian  current  and  wave 
propagation  direction  are  known. 

Spline  of  Record 

In  order  to  avoid  being  restricted  by  the  sampling  rate  of  a  particular  set  of 
measurements,  a  continuous  record  is  computed  by  spline  interpolation  among  the 
measured  points.  This  allows  the  window  widths  and  the  number  and  location  of 
samples  in  each  window  to  be  chosen  independently  of  the  sampling  frequency  of  the 
data.  Care  must  be  taken,  however,  to  be  sure  to  include  an  adequate  segment  of 
the  record  in  each  window.  A  window  that  includes  only  a  single  measured  point  is 

computationally  possible,  but  would  result  in  misleading  accuracy  in  the  resultant 
solution. 
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Output  Locations 

The  spacing  of  the  desired  output  locations  must  be  chosen  to  determine  the 
placement  of  the  windows.  The  solution  in  each  window  is  expected  to  be  most 
accurate  in  the  center,  and  thus  a  separate  window  is  centered  at  each  location  where 
output  is  desired.  As  each  window  is  computed  independently,  there  is  no  restriction 
on  the  spacing  of  the  output  windows. 

Non-Dimensionalization 

The  comparisons  of  the  errors  in  equations  of  different  dimensional  quantities  may 
result  in  spurious  solutions  that  over  emphasize  certain  equations  in  the  formulation. 
While  the  familiar  dimensional  forms  of  the  equations  have  been  presented  here  for 
clarity,  all  parameters  and  variables  are  non-dimensionalized  by  physically  identifiable 
parameters  before  computation.  The  mass  density  of  water  (p),  acceleration  of  gravity 
(p),  and  the  mean  zero  crossing  frequency  of  the  measured  record  (0;^)  are  used  to 
define  characteristic  length,  time  and  mass  scales. 

Length  scale  = 

Time  scale  =  l/w^  (3 

3 

Mass  scale  =  — ^ 

3.3.2  Optimization  Procedure 

The  primary  process  in  the  LFI  method  is  a  nonlinear  optimization  of  a  system 
typically  involving  up  to  14  algebraic  equations  in  12  unknown  parameters.  A  system 
as  complex  as  this  may  have  numerous  local  minima  that  result  in  spurious  solutions. 
The  best  way  to  avoid  these  solutions,  and  to  ensure  efficient  optimization,  is  to  start 
with  a  good  estimate  for  the  unknowns  in  the  system,  and,  in  addition,  to  have  a  set 
of  criteria  for  identifying  spurious  solutions. 
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Initial  Estimate 


The  first  step  in  each  window  is  to  establish  an  initial  estimate  for  the  optimization 
procedure.  Linear  wave  theory  can  be  used  to  produce  estimates  for  the  parameters 
of  correct  magnitudes.  The  linear  subsurface  dynamic  pressure  is: 


Pd  =  a  cos  {kx  —  u)t) 

cosh  kih  +  2p)  (^•^) 

cosh  kh 

where  A.  is  the  amplitude  of  the  linear  potential  function.  The  wave  number,  A;,  and 
frequency,  uj,  are  related  through  the  linear  dispersion  relation: 


(a;  —  kUy  =  gk  tanh  kh  (3.9) 

Frequency  Nielsen  (1986,  1989)  established  a  method  for  determining  the  param¬ 
eters  of  a  local  linear  approximation  to  waves  from  a  pressure  record.  His  method 
determined  a  local  frequency  and  wave  number  that  could  be  used  to  find  the  location 
of  the  water  surface.  A  similar  method  is  used  here  to  determine  the  first  estimate 
for  the  local  frequency  and  wave  number.  Frequency  of  a  sinusoidal  signal  of  the  form 
Pd  =  a  cos  (kx  —  ut)  is  available  from  the  second  derivative: 


I  1  d^pd 

This  approach  requires  an  estimate  for  the  value  of  the  second  time  derivative  of 
the  record  throughout  each  window.  An  estimate  for  the  second  derivative  is  directly 
available  from  the  cubic  spline  of  the  measured  points.  This  estimate,  however,  is 
very  sensitive  to  random  error  in  the  measurements.  Nielsen  suggests  estimating  the 
value  of  the  derivative  through  the  use  of  divided  differences.  This  solution  works 
well  for  a  smooth,  sinusoidal  record,  but  is  also  very  sensitive  to  noise  in  the  record, 
particularly  in  areas  of  small  curvature;  estimating  a  second  derivative  from  a  small 
segment  of  a  noisy  record  can  result  in  large  errors.  In  order  to  accommodate  the 
inevitable  noise  in  an  actual  record,  a  different  approach  is  taken  here. 


In  each  window,  a  third  order  polynomial  is  fit  to  the  record  in  the  least  squares 
sense.  The  second  derivative  throughout  the  window  can  then  be  computed  from  this 
polynomial.  By  using  more  than  four  points  in  each  window,  any  noise  in  the  record 
is  smoothed  out  by  the  least  squares  fit.  This  approach  proved  to  be  consistently 
reliable  for  artificially  generated  noisy  records,  resulting  in  reasonable  estimates  for 
the  value  of  the  second  derivative  throughout  the  window.  A  set  of  frequencies  is 
computed  from  the  estimate  of  the  second  derivative  at  each  of  the  considered  nodes. 
The  mean  of  these  frequencies  is  used  as  a  first  estimate  for  the  local  frequency  of  the 
record. 


Amplitude  and  Phase  Once  the  frequency  is  known,  the  amplitude  and  spatial 
phase  (kx)  of  a  particular  record  can  be  found  b}'  rearranging  the  equation  as  a  linear 
least  squares  problem  by  separating  the  cosine  and  sine  components: 


Pd.i  =  acos{kx  —  uti)  =  6cosu.’t,  +  csinwt,  (3-11) 

where  a  =  \/b'^  +  and  kx  =  arctan  (c/6).  This  results  in  the  following  matrix 

equation  in  the  unknown  amplitudes,  h  and  c. 

Pd,i 

Pd.2 

(3.12) 

PdJ. 

The  system  is  determined  if  pd{t)  is  defined  at  two  points.  If  /  >  2,  the  system  is 
over  determined  and  can  be  solved  in  the  least  squared  sense  by  common  algorithms 
in  numerical  linear  algebra  libraries.  This  method  provides  a  first  estimate  for  the 
parameters:  a,  kx,  and  oj. 


Refining  the  Linear  Estimates  These  estimates  for  the  parameters  can  be  quite 
poor,  as  they  are  all  dependent  on  the  initial  estimate  for  the  second  time  derivative 
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of  the  record.  In  order  to  improve  them,  the  estimates  are  refined  by  optimizing  for 
the  best  least  squares  fit  in  the  window. 


I 

minimizeO(X)  =  ^  -  a  cos  {kx  -  uU))^ 

i=l 

X  =  (a,  kx^io) 


(3.13) 


Where  is  the  measured  dynamic  pressure  at  time  U.  The  result  of  this  optimiza¬ 
tion  is  a  linear  estimate  for  the  pressure  record  that  fits  the  measured  record  most 
closely  in  the  window. 


Wave  Number  and  Fourier  Amplitudes  Once  the  optimum  initial  frequency, 
phase  and  amplitude  are  found,  the  wave  number  is  computed  with  the  linear  disper¬ 
sion  relation  (Eq.  3.9),  and  the  first  estimates  for  Fourier  amplitudes  are  computed 
as  follows: 


a  cosh  kh 
puj  cosh  A;(h  -|-  Zp) 


(3.14) 


Aj  =  qMi  (3.15) 

The  value  of  a  is  not  critical  and  a  =  0.1  was  found  to  be  satisfactory. 


Water  Surface  The  location  of  the  water  surface  at  the  N  nodes  throughout 
the  window  is  estimated  from  the  linear  pressure  response  function  with  stretching 
(Nielsen  1989): 


^  Pdiin)  cosh  {k  {h  -I-  Pd{tn)/P9)) 
pg  cosh  k{h  +  Zp) 


(3.16) 


where  Zp  is  the  vertical  location  of  the  pressure  gauge,  and  is  the  elevation  of  the 
water  surface  node  at  pd{t)  is  computed  from  Eq.  3.8. 
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Nonlinear  Optimization 

Once  there  is  a  reasonable  first  estimate  for  all  the  unknown  values  in  the  potential 
function,  standard  nonlinear  optimization  routines  are  adequate  for  this  system.  For 
the  results  given,  the  Levenberg-Marquardt  algorithm  was  used  as  implemented  by 
the  Matlab  Optimization  Toolbox  (Grace  1992). 

If  the  optimization  routine  successfully  finds  a  minimum,  the  solution  is  checked 
to  see  if  a  clearly  spurious  solution  is  found.  Spurious  solutions  are  identified  by  the 
following  criteria: 

•  Very  large  or  highly  variable  errors. 

•  First  order  amplitude  smaller  than  higher  order  amplitudes. 

•  Unrealistically  large  or  small  frequency  or  wave  number. 

•  Large  discontinuity  between  windows  in  the  predicted  water  surface. 

It  is  unusual  for  the  optimization  routine  to  converge  to  a  spurious  solution.  It  is  far 
more  common  for  the  routine  not  to  converge  at  all. 

If  no  solution  is  found,  or  a  spurious  solution  is  found,  it  is  necessary  to  revise  the 
parameters  of  the  optimization.  For  the  next  attempt,  the  window  width  is  increased 
by  a  factor  of  1.5  (1.5ro),  and  the  procedure  is  repeated.  If  this  is  not  successful, 
the  window  width  is  increased  once  more  to  twice  the  primary  width  (2ro).  When 
increasing  the  window  width  is  not  successful,  the  order  of  the  potential  function  is 
decreased  until  a  solution  is  found.  If  none  of  these  adjustments  result  in  an  acceptable 
solution,  the  window  is  skipped,  and  future  analysis  must  be  interpolated  through 
that  point.  These  adjustments  are  most  likely  to  be  needed  in  the  long,  flat  trough 
of  a  shallow  water  wave,  where  the  window  needs  to  be  expanded  to  include  some 
curvature  to  indicate  the  frequency. 

Difficulties  may  occur  in  addition  near  zero  crossings,  where  there  is  also  little 
curvature  in  the  record.  Here  the  effects  of  amplitude  and  frequency  may  not  be 
independent,  as: 

lim  a  sin  (ut)  =  au>t 

t-fO 


(3.17) 
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At  the  limit  near  the  zero  crossing,  a,  and  lo  have  the  same  effect,  and  the  optimization 
routine  can  not  distinguish  them.  Widening  the  window  to  include  more  of  the 
surrounding  record  avoids  this  difficulty,  and  is  generally  successful  in  this  situation 
as  well  as  in  long,  flat  troughs. 

Another  complication  can  be  a  record  that  is  exactly  symmetric  about  the  crest 
of  a  wave.  In  this  case,  the  equations  on  either  side  of  the  crest  are  not  independent. 
This  situation  is  unlikely  to  arise  in  a  field  record,  and  can  easily  be  accommodated 
by  using  an  asymmetric  distribution  of  points  in  that  window. 

3.4  Theoretical  Records 

To  avoid  complications  from  measurement  error  in  the  initial  testing  of  the  method, 
pressure  records  generated  by  Fourier  Steady  wave  theory  (Sobey  1989)  were  used. 
This  approach  also  has  the  advantage  of  providing  a  solution  with  the  complete  kine¬ 
matics,  to  compare  with  results  from  the  LFI  method.  Field  or  laboratory  data  that 
includes  a  full  set  of  measured  kinematics  are  not  available.  Fourier  theory  provides 
a  near-exact  solution  for  irrotational  steady  waves  and  can  be  applied  at  any  depth 
(Rienecker  and  Fenton  1981;  Sobey  1989).  The  use  of  analytically  generated  records 
also  allows  the  method  to  be  tested  under  a  large  range  of  conditions,  by  varying 
water  depth,  wave  period,  and  wave  height.  This  is  useful  for  establishing  criteria  for 
choosing  the  solution  parameters,  such  as  window  width,  number  of  nodes  in  each 
window,  and  order  of  the  solution. 

Shallow  Water  To  demonstrate  the  optimization  procedure  in  a  window,  figure  3.2 
summarizes  the  results  from  an  initial  estimate,  before  the  final  optimization.  This 
is  a  window  near  the  crest  of  a  steep,  shallow  water  wave  generated  by  18th  order 
Fourier  theory.  The  parameters  of  the  wave  are:  5m  water  depth,  3m  wave  height,  10s 
period,  and  zero  Eulerian  current,  with  the  pressure  record  measured  at  the  bottom. 
The  parameters  of  the  LFI  solution  are:  sixth  order  (^J  =  6),  and  window  width  two 
seconds  (tq  =  0.2r^)  ,  centered  0.5s  before  the  crest. 

The  top  plot  shows  the  measured  dynamic  pressure  as  given  by  the  Fourier  theory 
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generated  record,  and  the  values  computed  from  the  Bernoulli  equation  (Eq.  3.4), 
with  the  parameters  of  the  potential  function  generated  by  the  initial  estimate.  The 
second  plot  is  the  water  surface  as  predicted  by  Fourier  theory  and  the  elevations  of 
the  water  surface  nodes  generated  by  the  initial  estimate.  Note  that  the  location  of 
the  actual  surface  is  given  in  the  plot,  but  it  is  not  available  to  help  determine  the 
solution.  These  points  were  all  generated  by  the  method  outlined  in  the  previous 
section,  with  only  the  pressure  record  as  a  guide. 

The  third  plot  shows  the  non-dimensional  errors  in  the  objective  functions:  the 
Bernoulli  equation  and  the  free  surface  boundary  conditions  (Eqs.  3.2,  3.3,  and  3.4). 
These  are  the  errors  that  are  minimized  by  the  optimization  to  find  the  solution.  If 
the  solution  were  perfect,  the  error  in  all  equations  throughout  the  window  would  be 
zero. 

The  initial  estimates  for  the  dynamic  pressure  and  the  water  surface  have  order  of 
magnitude  and  trend  agreement.  The  errors  in  the  objective  functions  are  on  order 
of  .03  and  show  a  systematic  pattern,  particularly  in  the  Bernoulli  equation.  This 
pattern,  and  the  fact  that  the  sharp  crest  of  the  wave  has  not  been  matched  indicate 
that  a  better  solution  can  be  found. 

The  results  after  the  nonlinear  optimization  are  given  in  Fig.  3.3.  At  this  point  the 
prediction  for  the  dynamic  pressure  is  essentially  exact.  This  is  virtually  always  the 
case,  as  the  pressure  record  is  available,  and  the  solution  is  optimized  to  that  record. 
The  predictions  for  the  water  surface  are  also  extremely  close.  This  is  an  impressive 
achievement,  as  location  of  the  water  surface  was  found  only  by  minimizing  the  errors 
in  the  free  surface  boundary  conditions.  The  non-dimensional  errors  in  the  Bernoulli 
equation  and  free  surface  boundary  conditions  are  on  order  of  .001,  and  show  no  clear 
trend.  The  lack  of  trend  indicates  that  the  remaining  error  is  random,  and  a  good 
solution  has  been  found.  The  LFI  method  was  able  to  capture  accurately  the  crest 
of  a  steep  shallow  water  wave. 

Figure  3.4  shows  the  results  of  the  method  for  the  complete  shallow  water  steady 
wave.  The  LFI  method  finds  the  water  surface  and  the  kinematics  on  the  surface 
essentially  exactly.  While  these  results  show  the  complete  wave,  each  of  the  indicated 
points  is  in  the  center  of  a  separate  window,  and  each  window  was  computed  com- 
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Figure  3.5:  Evolution  of  the  solution  for  a  shallow  water  wave 

pletely  independently  of  the  other  windows.  In  this  case,  the  parameters  of  the  LFI 
solution  are:  primary  window  width  of  2s  (tq  =  0.27^),  with  a  sixth  order  potential 
function  ( J  =  6),  seven  samples  on  the  pressure  record,  and  seven  water  surface  nodes 
^  resulting  in  21  equations  in  16  unknowns  in  each  window.  The  window 
width  of  two  seconds  is  one  fifth  of  the  period  of  the  wave,  and  is  a  reasonable  length 
of  time  to  extend  the  locally  steady  approximation. 

Evolution  of  the  Solution  Fig.  3.5  shows  the  evolution  of  the  solution  as  the 
window  is  moved  along  the  wave.  The  top  figure  shows  the  non-dimensional  values  of 
the  local  fundamental  wave  number,  k,  and  the  local  fundamental  frequency,  w.  The 
importance  of  the  local  nature  of  the  solution  is  apparent  in  this  figure,  as  the  solution 
varies  substantially  from  window  to  window.  The  wave  number  and  frequency  both 
increase  to  maximum  magnitude  near  the  crest  (4  <  tu;^  <  7),  and  have  lower  values 
in  the  trough  region  (near  =  3  and  tu)^  =  9).  This  pattern  suggests  that  the 
kinematics  in  the  crest  region  are  similar  to  those  of  a  much  higher  frequency  lower 


40 


order  wave,  and  the  trough  kinematics  similar  to  a  lower  frequency  wave. 

Deep  Water  Figure  3.6  shows  the  results  of  the  LFI  method  for  an  entire  deep 
water  wave.  The  wave  was  generated  by  10th  order  Fourier  wave  theor}',  with  param¬ 
eters:  100m  water  depth,  10m  wave  height,  10s  period,  and  zero  Eulerian  current, 
with  the  pressure  record  measured  10m  below  the  mean  water  level.  The  parameters 
of  the  LFI  solution  are:  Primary  window  width  of  Is  (tq  =  O.lT,),  wdth  a  fourth  order 
potential  function  ( J  =  4),  five  samples  on  the  pressure  record,  and  five  w'ater  surface 
nodes(/  =  N  =  5),  resulting  in  15  equations  in  12  unknowns  in  each  window.  Once 
again,  the  LFI  solution  matches  the  actual  solution  essentially  exactly. 

3.4.1  Choosing  Parameters  of  Solution 

Unlike  steady  w^ave  theory  where  only  the  order  requires  prior  specification,  this 
LFI  application  requires  prior  specification  of  four  parameters: 

•  Order  of  the  solution  ( J) 

•  Window  width  (tq) 

•  Number  of  nodes  on  the  water  surface  (N) 

•  Number  of  nodes  on  the  pressure  record  (/) 

Experimentation  h«us  demonstrated  that  the  resulting  solutions  are  not  highly  sensi¬ 
tive  to  the  values  of  these  parameters;  however,  a  reasonable  solution  will  not  result 
if  these  parameters  are  not  selected  judiciously.  The  experience  acquired  while  devel¬ 
oping  the  LFI  method  on  analytically  generated  records  has  provided  some  guidelines 
to  help  select  appropriate  values.  Despite  these  guidelines,  indiv'^idual  judgment  must 
be  used  when  applying  the  LFI  method  to  a  particular  measured  record. 

Number  of  Nodes  on  the  Water  Surface  and  Pressure  Record  Mathemat¬ 
ically,  the  system  of  equations  is  specified  if  there  are  at  least  as  many  equations  as 
unknown  parameters.  In  this  case,  /  -f  >  4  -1-  J.  While  meeting  this  criterion  is 
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Figure  3.6:  LFI  predictions  ( J  =  4)  and  analytical  kinematics  at  the  predicted  water 
surface  for  deep  water  wave 
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the  minimum  mathematical  requirement,  there  are  other  factors  that  must  be  taken 
into  account  to  assure  a  reasonable  physical  solution.  The  local  water  surface  is  rep¬ 
resented  by  N  nodes  in  each  window.  Defining  the  surface  with  A’  =  7  -f  1  points 
achieves  a  representation  for  the  water  surface  of  the  same  order  as  the  potential 
function  being  used.  Care  must  be  taken  with  v'ery  low  order  solutions;  for  example, 
2  points  would  define  the  water  surface  to  first  order,  but  would  provide  no  indication 
of  the  curvature  of  the  surface.  As  a  rule  of  thumb,  at  least  three  points  should  be 
used,  regardless  of  order. 

Defining  the  water  surface  at  J  +  1  points  provides  2(J  +  1)  equations  (/^  and 

,  at  the  J-f  1  points).  To  keep  the  order  of  approximation  consistent,  the  Bernoulli 
equation  is  also  applied  at  J-f  1  points  within  each  window.  Using  J  +  1  water  surface 
and  pressure  record  nodes  overspecifies  the  solution  at  all  orders.  This  arrangement 
of  equations  proved  effective  for  all  the  examples  given  in  this  chapter  on  analyticall}'^ 
derived  records. 

While  that  approach  was  effective  on  these  few  examples,  it  is  important  that  the 
final  solution  matches  the  measured  record  closely.  It’s  possible  that  this  suggested 
arrangement  of  points  w'ould  allow  the  optimization  routine  to  be  biased  towards 
the  more  numerous  free  surface  boundary  conditions,  giving  a  solution  that  does  not 
match  the  pressure  record  w'ell.  The  Bernoulli  equation  could  be  applied  at  more 
points  on  the  pressure  record  than  the  number  of  nodes  defining  the  water  surface. 
Increasing  the  number  of  points  at  which  the  Bernoulli  equation  is  applied  will  shift 
the  emphasis  of  the  optimization  to  the  measured  record.  Additional  points  on  the 
measured  record  can  also  be  useful  for  accommodating  measurement  noise  that  may 
be  present  in  field  or  laboratory  records. 

Order  of  Solution  Similarly  to  high  order  steady  wave  theory,  the  order  chosen 
for  the  solution  is  influenced  by  a  number  of  factors  including  the  height  of  the  waves, 
the  depth  of  the  water,  and  the  accuracy  desired.  As  with  steady  w'ave  theory,  higher 
order  results  in  greater  accuracy  at  the  expense  of  computational  simplicity,  and  is 
necessary  for  larger  waves  and  for  shallow  water.  The  following  examples  will  help  to 
provide  guidelines  for  the  order  chosen. 


Figure  3.7:  Crest  of  a  shallow  water  wave  at  order  1 


Figure  3.8:  Crest  of  a  shallow  water  wave  at  order  2 


Figure  3.9:  Crest  of  a  shallow  water  wave  at  order  3 


Figure  3.10:  Crest  of  a  shallow  water  wave  at  order  4 


Figure  3.12:  Crest  of  a  shallow  water  wave  at  order  6 
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Figures  3.7  through  3.12  show  the  crest  of  the  same  shallow  water  wave  showm  in 
Fig.  3.4,  computed  to  orders  J  =  1  through  6.  At  first  order,  the  non-dimensional 
errors  in  the  Boundary  conditions  and  the  Bernoulli  equation  are  of  order  10~®.  These 
are  small  errors,  but  the  sharp  curvature  at  the  crest  has  not  been  captured.  As 
the  order  is  increased,  the  magnitude  of  the  errors  increases.  The  increase  in  the 
magnitude  of  the  errors  is  due  to  the  increase  in  the  number  of  equations  included. 
There  are  three  additional  equations  included  in  the  optimization  for  each  increase 
in  order  (the  two  FSBCs  at  each  additional  water  surface  node,  and  one  additional 
pressure  record  node),  but  there  is  only  one  more  parameter  in  the  solution  vector. 
In  order  to  best  satisfy  this  system,  the  optimization  finds  a  solution  that  results  in 
slightly  more  error  in  each  equation.  The  magnitude  of  the  errors  does  increase,  but 
the  solution  slowdy  converges  to  very  precisely  match  the  sharp  crest  at  sixth  order. 
An  asymmetric  distribution  of  points  was  used  in  this  window  to  accommodate  the 
symmetry  about  the  crest. 

Figures  3.13  through  3.16  show  the  crest  of  the  same  deep  water  wave  as  Fig.  3.6, 
computed  to  orders  1  through  4.  Even  at  first  order.  Fig.  3.13,  the  solution  is  very 
good.  Deep  water  waves  generally  do  not  require  a  very  high  order  solution,  linear 
wave  theory  often  being  reasonably  adequate  in  these  conditions.  It’s  important  to 
keep  in  mind  that,  although  this  solution  is  first  order,  it  is  still  nonlinear,  having 
found  a  minimum  in  the  errors  of  the  full  nonlinear  governing  equations,  and  the 
frequency  and  w'ave  number  are  free  to  vary,  not  being  bound  by  the  linear  dispersion 
relationship.  The  local  nature  of  the  solution  would  allow  it  to  change  with  time, 
accommodating  an  irregular  profile  at  low  order  better  than  an  equally  low  order 
global  solution. 

In  this  case,  first  order  provides  an  acceptable  solution;  however,  the  water  surface 
is  more  accurately  matched  as  the  order  is  increased,  and  the  solution  converged  well 
at  higher  order,  so  there  is  little  penalty  in  using  up  to  fourth  order.  For  an  irregular 
record,  higher  order  is  more  likelj'  to  be  successful  in  matching  the  irregularity  in  the 
record.  As  the  examples  show,  deep  water  waves  can  be  well  represented  at  low  order, 
while  higher  order  in  necessary  in  shallow  water. 


Figure  3.13:  Crest  of  a  deep  water  wave  at  order  1 


Figure  3.14:  Crest  of  a  deep  water  wave  at  order  2 
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Window  Width  Choice  of  window  width  is  a  balance  between  including  sufficient 
curvature  in  the  record  to  identify  the  local  wave  response  while  minimizing  the  extent 
of  the  locally  steady  approximation.  Another  factor  to  be  considered  is  the  sampling 
rate  of  the  data.  The  LFI  method  employs  a  cubic  spline  for  interpolation  among 
the  measured  points.  This  approach  allows  the  selection  window  widths  and  sampling 
locations  without  being  restricted  by  the  data  sampling  locations.  While  this  freedom 
is  very  useful,  it  can  be  deceiving,  as  the  data  are  not  truly  continuous.  As  the  data 
are  not  continuous,  it  is  important  to  be  sure  that  the  window  width  chosen  includes 
sufficient  measured  data  points  to  justify  the  order  being  used.  Unfortunately,  field 
data  are  often  sampled  at  a  frequency  of  as  low  as  IHz.  In  this  case,  a  window  width 
of  one  tenth  the  mean  zero  crossing  might  be  on  order  of  1  second.  This  would  provide 
a  maximum  of  two  actual  data  points  in  each  window.  It  may  not  be  appropriate 
to  attempt  a  high  order  solution  with  such  little  hard  data.  Ideally,  the  LFI  method 
should  be  used  with  data  sampled  at  a  higher  rate.  If  that  isn’t  possible,  wider 
windows,  and  perhaps  lower  order  solutions,  should  be  used. 

In  the  case  of  the  theoretical  records  used  in  the  previous  section,  the  data  sam¬ 
pling  rate  could  be  arbitrarily  high.  Without  being  restricted  by  measurement  fre¬ 
quency ,  some  guidelines  for  window  width  have  been  determined.  A  primary  window 
width  of  one  tenth  the  zero  crossing  period  provides  a  good  starting  point.  Using 
a  smaller  window  often  did  not  allow  convergence  of  the  solution.  It  was  necessary 
to  occasionally  increase  window  width  at  some  locations  on  the  record,  as  discussed 
m  section  3.3.2.  When  a  primary  window  width  is  determined  that  works  for  most 
of  the  record,  while  needing  to  be  increased  at  a  only  a  few  locations,  the  optimum 
width  has  been  found. 

In  the  case  of  the  shallow  water  wave  discussed  in  the  previous  section,  a  window 
width  of  O.lT^j  was  successful  for  up  to  a  fourth  order  solution.  Unfortunately,  the 
fourth  order  solution  did  not  fully  capture  the  sharp  crest  of  the  wave.  When  a  higher 
order  solution  was  attempted,  the  optimization  would  not  converge  with  the  given 
window  width,  A  wider  window  had  to  be  used  in  order  to  obtain  a  solution  of  high 
enough  order  to  capture  the  sharp  crest  of  the  wave,  and  thus  a  window  width  of0.2Tz 
was  used  in  that  example.  In  general,  it  has  been  necessary  to  use  wider  windows  for 
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higher  order  solutions. 


3.5  Laboratory  Measurements 

The  following  results  use  laboratory  data  collected  by  Murray  Townsend  and 
John  Fenton  at  the  Australian  Maritime  Engineering  Cooperative  Research  Center  at 
Monash  University,  Melbourne,  Australia  in  June,  1996.  The  dynamic  pressure,  fluid 
velocities  and  water  surface  were  measured  at  the  same  horizontal  location  along  the 
flume;  the  pressure  at  a  variety  of  elevations,  and  the  horizontal  and  vertical  velocities 
at  an  elevation  of  -0.8  meters,  with  the  origin  at  the  mean  water  level.  The  still  water 
level  was  1.55  meters,  with  the  data  sampled  at  60  Hz. 

The  wave  flume  is  52m  long  by  2.2m  wide  with  two  working  sections  4m  and  2.5m 
long  connected  by  a  ramp.  A  false  floor  had  been  built  into  the  shallower  working 
section  to  bring  the  depth  to  1.55m.  At  one  end  of  the  flume  is  a  hydraulically  oper¬ 
ated  bottom  hinged  wave  paddle  and  at  the  other  a  beach  that  absorbs  a  minimum 
of  96%  of  wave  energy  across  the  range  of  operating  frequencies.  The  measurements 
were  taken  approximately  22m  from  the  beach  end  of  the  flume.  The  beach  length  is 
6m. 

The  wave  surface  was  measured  with  a  capacitance  wave  probe  with  typical  cali¬ 
brations  providing  a  correlation  coefficient  P?  greater  than  0.999.  The  pressure  was 
measured  with  Druck  PDCR35/D  transducers  located  in  the  flume  near  the  center 
of  a  long  (2.5m  long  by  2m  high)  plywood  board  aligned  with  the  direction  of  wave 
propagation.  The  measurement  face  of  the  probe  was  flush  with  the  board  surface 
to  eliminate  local  flow  separation.  The  values  from  calibrations  were  in  all  cases 
above  0.999.  Fluid  velocities  were  measured  with  a  factory  calibrated  Alec  electronics 
ACM250-A  electromagnetic  current  meter  which  is  capable  of  velocity  measurements 
in  two  dimensions. 

One  of  the  difficulties  encountered  with  pressure  measurements  is  the  question 
of  what  has  actually  been  measured.  If  the  transducer  itself  is  absolutely  accurate, 
the  pressure  recorded  wdll  be  the  actual  pressure  at  the  location  of  the  transducer. 
Unfortunately,  the  presence  of  the  instrument  is  likely  to  alter  the  flow  in  its  vicinity. 
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creating  substantial  differences  between  the  pressure  at  the  transducer,  and  the  pres¬ 
sure  that  would  exist  if  the  transducer  were  not  there  to  disturb  the  flow.  Another 
source  of  error  is  the  frequency  response  of  the  transducer  (Raichlen  et  al.  1990; 
Ippen  and  Raichlen  1957).  The  frequency  response  of  pressure  transducers  is  not  flat, 
and  as  a  result,  the  recorded  signal  could  be  quite  different  from  the  actual  pressure. 

In  the  experimental  setup  described  above,  the  pressure  transducers  were  mounted 
on  a  large  plywood  panel  oriented  in  line  with  the  flume.  The  measurement  face  of 
the  instrument  was  flush  with  the  surface  of  the  plywood  to  minimize  dynamic  effects 
near  the  gauge.  The  panel  was  large  enough  for  the  boundary  layer  to  be  fully 
developed  near  the  surface  of  the  plywood,  to  prevent  edge  effects  from  the  edge 
of  the  plywood  affecting  the  measurements.  This  arrangement  is  expected  to  have 
resulted  in  minimal  dynamic  effects  on  the  measured  pressure. 

There  is  no  information  available  about  the  frequency  response  of  the  pressure 
transducers  used  for  this  experiment.  However,  a  spectral  analysis  of  the  records  can 
help  identify  any  potential  problems.  The  top  plot  of  figure  3.17  gives  the  measured 
water  surface  and  subsurface  pressure  for  a  short  segment  of  a  record.  There  are 
some  clear  higher  frequency  fluctuations  in  the  water  surface  that  do  not  appear  in 
the  pressure  record.  In  this  particular  segment,  there  is  a  sharp  secondary  crest  in 
the  first  trough  (near  the  2s  point)  in  the  water  surface  record. 

This  loss  of  high  frequency  information  in  the  pressure  records  is  confirmed  by 
an  examination  of  the  discrete  Fourier  transform  of  the  records.  The  second  plot  in 
figure  3.17  is  the  smoothed  variance  spectra  of  the  two  records  from  which  the  above 
segments  were  taken.  The  water  surface  record  has  a  great  deal  more  variance  at 
the  higher  frequencies.  Note  that  the  spectrum  of  the  pressure  record  has  decayed  to 
almost  zero  at  a?  =  4s  ^  point,  while  there  is  still  considerable  variance  in  the  water 
surface  at  that  frequency. 

Linear  wave  theory  predicts  that  the  motion  of  high  frequency  waves  decays  with 
depth  much  faster  than  low  frequencies.  In  this  example,  the  non-dimensional  water 
depth  at  the  peak  frequency  is  u'^h/g  fs  1.  This  is  generally  considered  intermediate 
depth  water,  and  the  decay  with  depth  is  expected  to  be  moderate.  At  the  frequency 
where  there  is  essentially  no  energy  in  the  pressure  record,  the  non-dimensional  water 
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depth  is  iJ^hlg  fa  2.5,  which  is  considered  to  be  deep  water.  In  deep  water,  the 
action  of  the  waves  decays  very  quickly  with  depth,  so  that  little  energy  is  expected 

to  remain  at  half  the  water  depth,  the  depth  at  which  the  pressure  measurements 
were  taken. 


The  third  plot  in  figure  3.17  is  the  spectrum  of  the  subsurface  pressure  record. 
The  solid  line  was  computed  from  the  measured  record,  and  the  dashed  line  was 
computed  from  the  water  surface  record,  using  the  linear  pressure  response  factor: 


cosh{kh) 


(3.18) 


where  ap{u;)  and  a„(w)  are  the  magnitude  of  the  pressure  and  water  surface  amplitude 
spectra  at  a  given  frequency,  and  the  wave  number,  k,  is  computed  from  the  linear 
dispersion  relationship: 


{uj  —  kUy  =  gktmh(kh)  (3.19) 

The  predicted  spectrum  matches  the  measured  spectrum  very  closely.  This  indicates 
that  the  loss  of  high  frequency  information  in  the  subsurface  pressure  record  is  the 
result  of  the  predicted  decay  with  depth  of  the  energy  at  the  higher  frequencies,  rather 
than  a  result  of  the  frequency  response  of  the  pressure  gauge  used. 

The  final  plot  in  figure  3.17  is  the  spectrum  of  the  water  surface  record.  The 
solid  line  was  computed  from  the  measured  water  surface,  and  the  dashed  line  was 
computed  from  the  pressure  record,  also  using  the  linear  pressure  response  factor.  The 
predicted  spectrum  of  the  water  surface  matches  the  measured  spectrum  well  near 
the  peak  frequency,  but  strongly  over-predicts  the  high  frequencies.  If  this  method 
were  used  to  predict  the  water  surface,  a  somewhat  arbitrary  high  frequency  cut  off 
would  have  to  be  established. 

An  examination  of  the  top  plot  in  figure  3.17  reveals  what  appears  to  be  two 
dominant  free  modes  in  the  water  surface.  There  is  a  large  low  frequency  mode 
with  a  period  of  approximately  2s,  and  a  higher  frequency  mode  with  a  period  of 
approximately  Is.  If  the  LFI  method  were  applied  to  the  water  surface  record  (Sobey 
1992),  the  moving  widow  could  capture  a  different  dominant  mode  in  each  window. 
Unfortunately,  the  higher  frequency  mode  has  decayed  with  depth  too  much  to  be 
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Water  Surface 


Figure  3.18:  Measured  and  predicted  water  surface  near  the  crest  of  a  sharp  wave  in 
the  flume 

clear  in  the  pressure  record;  applying  the  LFI  method  to  the  pressure  record  results 
in  predictions  that  capture  the  dominant  low  frequency  mode,  and  miss  the  higher 
frequency  mode.  If  the  higher  frequency  information  is  not  in  the  local  segment  of 
the  record,  it  is  not  possible  to  reconstruct  it. 

In  the  previous  section,  the  LFI  method  was  applied  to  pressure  records  generated 
by  Fourier  steady  wave  theory.  These  records  had  only  a  single  free  mode,  and  thus 
the  LFI  method  was  able  to  identify  that  mode.  In  addition,  in  shallow  water,  there 
is  little  decay  with  depth  of  the  wave  action,  and  the  LFI  method  could  be  expected 
to  be  effective,  as  it  was  on  the  Fourier  record.  In  deep  water,  the  method  can  be 
effective  if  the  pressure  gauges  are  located  near  enough  to  the  surface.  For  the  deep 
water  steady  wave  presented  previously  (figure  3.6)  the  pressure  record  was  taken 
from  a  depth  of  one  wave  height  below  the  water  surface,  one  tenth  of  the  total 
depth.  At  this  depth,  there  is  adequate  information  to  identify  the  action  at  a  wide 
range  of  frequencies. 

Given  the  limitations  of  the  data,  a  segment  of  the  record  was  chosen  in  which 
there  did  not  appear  to  be  a  substantial  high  frequency  component  to  the  water 
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Figure  3.19:  Segment  of  record  from  flume  measurements 
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surface  record,  and  thus  the  effect  of  the  high  frequency  decay  should  cause  few 
problems,  and  the  LFI  method  can  be  applied. 

While  this  particular  segment  was  chosen  to  minimize  the  effects  of  the  high 
frequency  decay,  it  was  not  possible  to  eliminate  it  entirely.  Figure  3.18  shows  the 
measured  water  surface  and  the  LFI  predictions  of  the  water  surface  in  three  windows 
in  the  vicinity  of  a  sharp  crest  in  the  record.  The  solution  parameters  are:  J  =  3, 
To  =  0.3Tz,  and  I  =  N  =  4.  Both  the  pressure  and  the  velocities  were  measured 
at  an  elevation  of  0.8m  below  the  MWL,  w'ith  the  mean  water  depth  1.55m.  The 
LFI  method  located  the  water  surface  fairly  accurately,  although  it  clearly  missed 
the  very  sharp  peak  of  the  crest.  Varying  the  solution  parameters  did  not  improve 
the  accuracy  at  the  crest.  In  fact,  a  smaller  w’indow  width  resulted  in  a  discontinuity 
in  the  water  surface  between  w’indows.  The  inability  of  the  LFI  method  to  capture 
the  sharp  crest  is  likely  due  to  the  missing  high  frequency  information.  If  the  high 
frequency  part  of  the  signal  is  missing  from  the  examined  pressure  record,  there  is  no 
way  to  recapture  that  information.  In  order  for  the  method  to  be  more  effective  in 
this  situation,  the  pressure  would  have  to  be  measured  at  a  depth  closer  to  the  water 
surface. 

Figure  3.19  show's  the  results  of  the  LFI  method  applied  to  the  entire  segment  of 
the  record  from  which  the  crest  was  taken.  The  w-ater  surface  is  captured  fairly  well, 
as  are  both  the  vertical  and  horizontal  velocities. 


3.6  Discussion 

The  given  results  demonstrate  the  potential  of  the  LFI  method  in  the  interpre¬ 
tation  of  submerged  pressure  traces  in  a  variety  of  conditions.  In  the  case  of  steady 
waves,  the  LFI  method  accurately  computed  the  detail  of  the  w’ave,  using  only  data 
from  a  small  w'indow’  in  time.  In  particular,  the  method  w'as  able  to  capture  the 
pronounced  sharp  crest  of  a  steep,  shallow  water  wrave. 

The  method  did  not  perform  as  w'ell  on  laboratory  records,  failing  to  capture 
some  of  high  frequency  detail  in  the  w'ater  surface.  This  is  due  to  the  limitations  of 
subsurface  pressure  records,  w'here  much  of  the  high  frequency  information  is  missing 


57 


from  the  record.  The  method  is  likely  to  perform  better  in  shallow  water,  or  with 
records  that  are  measured  closer  to  the  water  surface.  Despite  this  limitation,  the 
method  was  able  to  capture  much  of  the  detail  of  a  irregular  laboratory  record,  and 
provide  fairly  accurate  estimations  of  the  local  kinematics. 

The  analysis  of  regular  waves  provides  guidelines  for  the  parameters  to  be  used 
in  the  analysis  of  irregular  waves.  In  shallow  water,  higher  order  solutions  and  wider 
windows  must  be  used  than  in  deep  water.  Window  widths  of  one  fifth  of  the  zero 
crossing  period  and  a  sixth  order  potential  function  are  adequate  for  the  shallowest 
waves,  and  window  widths  as  small  as  one  tenth  of  the  zero  crossing  period  and  a 
third  order  potential  function  are  adequate  for  deep  water. 

The  laboratory  results  indicate  that  the  method  requires  good  precision  and  care  in 
the  measurements.  Any  high  order  method  demands  very  accurate  data,  but  this  need 
IS  exacerbated  by  the  ill  posed  problem  of  determining  the  near  surface  kinematics 
from  a  subsurface  record.  W^hile  not  particularly  sensitive  to  random  noise  in  the 
record,  the  decay  with  depth  of  the  high  frequency  information  makes  it  very  difficult 
to  capture  the  high  frequency  modes.  Fundamentally,  the  LFI  method  is  designed  to 
capture  the  dominant  free  mode  in  each  given  window.  As  the  higher  frequency  modes 
decay  faster  with  depth,  and  if  the  measurements  are  taken  far  below  the  surface,  the 
dominant  mode  will  always  be  one  of  lower  frequency  modes.  This  difficulty  would 
be  exacerbated  by  any  limitations  in  the  frequency  response  of  the  gauges. 

The  computer  resources  required  for  the  method  are  substantial,  but  not  pro¬ 
hibitive.  As  computers  continue  to  get  faster,  computation  time  is  not  the  considera¬ 
tion  that  it  once  was.  The  method  was  developed  and  all  computation  done  using  the 
Matlab  computational  package.  Matlab  is  an  interpreted  language  that  provides 
a  very  easy  to  use  interface  to  a  robust  and  complete  library  of  computational  and 
visualization  routines,  allowing  for  rapid  development  of  new  methods.  Being  an  in¬ 
terpreted  language,  the  resulting  code  is  not  as  fast  as  it  might  be  if  the  routines  were 
programmed  in  a  compiled  language,  such  as  Fortran  or  C.  The  computational  speed 
is  also  affected  by  the  degree  to  which  the  program  pauses  to  provide  visual  output. 
Perhaps  a  better  measure  of  the  computational  intensity  of  the  method  is  a  count  of 
the  total  number  floating  point  operations  (flops)  used  to  compute  the  solution. 
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The  examples  in  this  dissertation  were  all  computed  on  a  Intel  90MHz  Pentium 
processor  based  PC  with  32MB  of  RAM,  running  the  Linux  operating  system.  As  an 
order  of  magnitude  estimate,  the  shortest  computation  time  for  an  example  in  this 
chapter  is  4.8  minutes  and  7.8  x  10®  flops  for  the  complete  laboratory  record  presented 
in  figure  3.19.  The  longest  computation  time  required  was  71  minutes  and  582  x  10® 
flops  for  the  shallow  water  wave  presented  in  figure  3.4.  There  are  a  number  of 
reasons  for  this  large  variation  in  computation  time.  The  first  is  simply  the  number 
of  window  solutions  computed.  As  each  window  solution  is  computed  separately, 
the  computation  time  increases  linearly  with  the  number  of  windows  computed.  If 
computational  time  is  a  concern,  this  can  be  taken  into  account  when  choosing  the 
output  spacing.  The  other  reason  that  the  shallow  water  wave  takes  much  longer 
to  compute  is  that  there  are  a  number  of  windows  that  did  not  converge  with  the 
first  set  of  computational  parameters.  The  optimization  is  run  for  many  iterations 
to  ensure  that  it  won’t  converge.  As  the  parameters  are  varied,  the  computation  is 
repeated.  This  process  takes  a  great  deal  of  time. 

With  further  development,  it  may  be  possible  to  determine  a  set  of  criteria  for  the 
computational  parameters  by  examining  the  segment  of  the  record  in  a  given  window. 
This  would  be  far  more  efficient  than  simply  attempting  a  solution  with  a  variety  of 
values  until  convergence  is  achieved. 
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Chapter  4 

Three  Dimensional  LFI  Theory 

The  previous  chapters  were  concerned  with  the  determination  of  wave  kinematics 
from  the  analysis  of  a  time  series  of  a  single  physical  quantity  at  a  single  location.  By 
assuming  that  the  flow  is  two  dimensional,  a  reasonable  approximation  of  the  wave 
kinematics  can  be  found.  Unfortunately,  this  analysis  does  not  give  any  information 
about  the  directionality  of  the  sea  state.  There  axe  some  processes  in  which  the  wave 
directions  are  directly  important,  such  as. sediment  transport.  Even  in  situations 
where  the  wave  directions  are  not  directly  important,  it  has  been  shown  that  omitting 
the  directional  nature  of  the  sea  results  in  substantial  errors  in  the  prediction  of  scalar 
quantities,  such  as  the  maximum  velocities  and  accelerations  in  a  measured  wave  crest, 
or  the  resultant  forcing  on  structures  (Dean  1977;  Forristall  et  al.  1978). 

In  order  to  capture  the  directional  nature  of  the  sea,  an  array  of  instruments  must 
be  used.  The  result  is  a  set  of  time  series  of  a  single  physical  quantity  at  a  number 
of  different  locations,  or  a  set  of  different  physical  quantities  at  the  same  or  different 
locations.  This  chapter  outlines  a  method  for  determining  the  directional  kinematics 
of  irregular  seas  that  can  be  adapted  to  accommodate  virtually  any  combination  of 
such  time  series. 
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seabed 

Figure  4.1:  Coordinate  system  used  for  3-d  method 


4.1  Three  Dimensional  Seas 

The  development  of  the  two  dimensional  LFI  method  was  closely  anchored  to 
the  very  complete  understanding  of  two  dimensional  steady  waves.  In  contrast,  the 
understanding  of  three  dimensional  wave  fields  is  not  nearly  as  complete.  Much  of 
the  literature  on  three  dimensional  seas  attempts  to  describe  the  motion  through  the 
use  of  a  directional  energy  spectrum.  Far  less  attention  has  been  directed  to  the 
determination  of  the  detailed  kinematics  of  directional  seas. 


4.2  Problem  Formulation 


The  governing  equations  for  three  dimensional  gravity  waves  are  a  straightforward 
extension  of  those  in  two  dimensions  to  include  the  third  dimension.  The  flow  is  taken 
to  be  irrotational  and  incompressible,  and  thus  the  kinematics  can  be  represented  by 
a  potential  function,  ^(x,j/,  2,t),  in  a  Cartesian  coordinate  system  (Fig.  4.1),  where: 


dx  ^  dy 


(4.1) 


u  and  V  are  the  horizontal  velocities  in  the  x  and  y  directions,  and  w  is  the  vertical 
velocity. 
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Field  Equation  The  field  equation  is  mass  conservation  for  irrotational  flow,  rep¬ 
resented  by  the  Laplace  equation: 


vV 


d^4> 

dx^  dy^  dz^ 


(4.2) 


Boundary  Conditions  The  boundary  conditions  are  the  bottom  boundary  condi¬ 
tion  (BBC)  for  a  horizontal  bed, 

"^  =  ^  =  0  at  z  =  -h  (4.3) 

the  dynamic  free  surface  boundary  condition  (DFSBC), 

.  1/  2  2  2x  - 

+  V  -f  w^)  +  grf  —  B  =  0  at  z  —  ij  (4.4) 

and  the  kinematic  free  surface  boundary  condition  (KFSBC), 
dr}  dr}  dr} 

¥  +  +  =  "  =  •'  (4-5) 

where  rf  is  the  elevation  of  the  water  surface. 

As  with  the  two  dimensional  formulation,  the  kinematic  free  surface  boundary 
condition  requires  the  gradients  of  the  water  surface  (|f,  g,  and  g).  Knowledge 
of  these  gradients  is  likely  to  be  limited.  In  the  case  of  an  array  of  water  surface 
measurements,  estimates  for  the  gradients  could  be  computed  by  interpolating  among 
the  measured  elevations,  but  this  would  result  in  a  low  order  estimate,  and  compound 
the  inevitable  error  in  the  measurements.  In  the  case  of  subsurface  measurements, 
there  are  no  data  as  to  the  location  or  the  gradients  of  the  water  surface. 

To  eliminate  these  gradients,  a  modified  kinematic  free  surface  boundary  condition 
is  defined  by  differentiating  the  DFSBC  following  the  motion  (Longuet-Higgins  1962), 
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as  was  done  in  two  dimensions  in  Section  2.1. 

MKFSBC  =  -p(KFSBC)  +  ^  (DFSBC)  = 


d^<f> 

dP 


av  av  du\ 

2^1/  dv  du' 

+  «  ^  +  uv—  +  uw— 
ox  ax  ax 

du  2  dw 

+  uv—  +  t»  -5-  +  vw— 
dy  dy  dy 

du  dv  2 
+  uw—  +  wv—  +  ic  =  0 


at 


z  =  r] 


(4.6) 


This  form  does  not  include  the  gradients  of  the  water  surface,  and  all  the  terms  are 
defined  by  the  potential  function.  Applying  both  the  modified  kinematic  free  surface 
boundary  condition  and  the  dynamic  free  surface  boundary  condition  completes  the 
formulation  without  the  need  for  knowledge  of  the  gradients  of  the  water  surface. 


Observational  Equations  The  field  equation  and  the  boundary  conditions  de¬ 
scribe  a  free  surface  potential  flow.  In  stead}'  two  dimensional  wav'e  theory,  a  wave 
height  and  periodicity  in  space  and  time  are  specified  to  complete  the  formulation. 
In  the  two  dimensional  LFI  theory  (Chapter  2),  a  form  for  the  potential  function  is 
chosen,  with  the  parameters  determined  to  fit  the  free  surface  boundary  conditions 
and  the  measured  record  in  a  small  window  in  time.  In  three  dimensions,  the  pro¬ 
cess  is  essentially  the  same.  A  three  dimensional  form  for  the  potential  function  will 
be  presented,  with  parameters  that  are  found  to  fit  the  measured  records  and  the 
boundary  conditions. 

In  order  to  define  a  solution  that  fits  the  measured  record,  observational  equations 
are  established.  These  equations  are  defined  to  make  use  of  the  particular  quantities 
that  have  been  measured.  In  the  case  of  an  array  of  water  surface  measurements,  they 
are  the  free  surface  boundary  conditions,  applied  at  the  measured  elev'ations  and  hor¬ 
izontal  locations  at  a  number  of  points  in  time  throughout  the  window  (Chapter  5). 
In  the  case  of  an  array  of  pressure  measurements,  they  are  the  Bernoulli  equation, 
applied  at  the  ele\'ation  and  horizontal  locations  of  the  measurements,  and  also  at 
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a  number  of  points  in  time  within  the  window  (Chapter  6).  The  method  could  be 
adapted  to  virtually  any  combination  of  measured  physical  quantities  by  establishing 
appropriate  observational  equations.  In  addition  to  the  aforementioned  water  surface 
and  pressure  measurements,  these  could  include  water  surface  gradient  or  accelera¬ 
tions  measured  by  wave  buoys,  subsurface  velocity  measurements,  or  any  combination 
of  these. 


4.3  Formulation  of  the  Solution  in  each  Window 

4.3.1  Background 

Two  dimensional  steady  wave  theory  provided  the  inspiration  for  the  development 
of  the  two  dimensional  LFI  method  (Chapter  2).  Unfortunately,  the  literature  does 
not  provide  as  solid  a  basis  for  nonlinear  interpretations  of  directional  seas.  There 
has,  however,  been  some  work  that  can  be  used  as  a  basis  for  a  directional  LFI 
method.  In  an  early  attempt  to  explore  the  nonlinear  nature  of  directional  seas, 
Longuet-Higgins  (1962)  computed  the  interaction  of  two  intersecting  steady  waves 
in  deep  water  through  the  use  of  a  double  perturbation  expansion  in  the  steepness 
of  the  waves  up  to  third  order.  The  result  was  a  potential  function  that  contained 

terms  representing  the  higher  order  interaction  between  the  phases  of  the  intersecting 
waves: 

4'^^^  —  Si  ^2)  -f  Si  —  S2) 

=  /.(Z.  25i  +  52)  +  /6(^.  2S.  -  S2)  +  /,(z,  S.  +  2S2)  +  /s(z,  5,  -  252) 

'S'n  =  (k„  •  X  -  Unt  Gn) 

where  n  [1,2],  1  is  the  mth  order  potential  function.  Si  and  S2  are  the  phase 

functions  of  the  two  waves,  is  the  vector  wave  number  of  wave  n,  x  is  the  horizontal 
position  vector,  and  a„  are  the  angular  frequency  and  initial  phase  of  the  nth 
wave.  The  functions,  /,,  were  determined  algebraically  by  expanding  the  free  surface 
boundary  conditions  in  a  Taylor  series  about  the  mean  water  level,  and  solving  for 


the  potential  function  at  each  order.  Nonlinear  frequency  modulation  was  not  taken 
into  account,  so  the  frequencies  and  wave  numbers  of  each  wave  independently  satisfy 
the  linear  dispersion  relation  for  water  of  infinite  depth: 


=^|kn| 


and  /2  are  independent  of  each  other,  and  are  the  familiar  linear  wave  solution: 


/i  =  sin  5i  sin  S2 

ki  /C2 


The  higher  order  terms  are  all  functions  of  both  waves,  and  thus  include  the  interac¬ 
tion  of  the  two  waves. 


A  number  of  investigators  subsequently  expanded  upon  this  work.  Hsu  and  Chen 
(1992)  presented  a  detailed  examination  of  Longuet-Higgins  (1962),  pointing  out  diffi¬ 
culties  that  arise  from  the  assumption  of  the  linear  dispersion  relation.  They  presented 
a  more  mature  analysis,  including  higher  order  modulation  of  the  wave  frequencies, 
and  higher  order  self  interactions  of  the  individual  waves.  This  resulted  in  a  com¬ 
plete  theory  up  to  third  order  for  two  intersecting  waves  in  deep  water.  Hsu  and 
Chen  also  proposed  a  systematic  ordering  in  the  phase  relationships  generated  by  the 
interactions  of  two  steady  waves  to  arbitrary  order. 

Expanding  upon  this  work,  Ohyama,  Jeng  and  Hsu  (1995a)  extended  the  per¬ 
turbation  expansion  method  in  a  number  of  ways.  The  most  recent  version  of  the 
method  accommodates  any  number  of  waves,  allows  for  water  of  finite  depth,  and  is 
accurate  to  fourth  order.  This  last  method  can  compute  the  water  surface  and  full 
kinematics  of  a  highly  irregular  sea,  as  produced  by  a  large  number  of  intersecting 
waves.  A  more  detailed  discussion  of  the  method  is  given  in  Appendix  A. 

Ohyama,  Jeng  and  Hsu’s  work  suggests  a  form  for  the  potential  function  that  can 
be  used  for  a  local  method  for  the  recreation  of  three  dimensional  kinematics  from 
arrays  of  wave  measurement  devices,  including  water  surface  arrays,  pressure  arrays, 
directional  current  meters,  or  any  combination  of  these. 

The  general  form  for  the  potential  function  representing  N  intersecting  steady 
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waves  at  order  J  is: 


^{x,  y,  z,  t)  =  UxX  +  Uyy  + 
j 


cosnAn(n  +  ^) 

cosh^llfe  "‘”'"11 


J-lill 

E  E 


[]  —  (*l5*2,-‘-  5*iv) 


<T[]  —  {iiSi  +  i2S2  +  •  •  •  +  iNSj^) 


Sn  {kx,n^  "t"  ky^fiy  Uint  -)- 


dy 


where  Ux  and  Uy  are  the  components  of  the  depth  uniform  Eulerian  current,  kx,n  and 
ky^n  are  the  vector  wave  number  components  of  the  nth  wave  in  the  x  and  y  directions. 
There  is  a  separate  summation  for  each  wave  considered. 

As  the  notation  in  Eq.  4.10  is  quite  confusing,  examples  with  two  and  three  waves 
are  given  below.  With  two  waves  {N  =  2): 


(f>{x,  y,  z,  t)  =  UxX  +  Uyy  + 


E 


2l  =  — «7 


■  ^-Nil 

,  E 

«2=“*/+|n  I 


tl,t2 


cosh  +  z) 

cosh  Ki^^i^h 


sm{iiSi  +  i2S2)  (4.11) 


Sn  {kx.n^  “f  ky_ny  ^nt  H"  ^n)  ^ 

~  yj {nkx^i^  +  *2^2;, 12)^  +  +  *2^2/, 


And  with  three  waves  {N  =  3): 


4>{x,  y,  2,  t)  =  UxX  +  Uyy  + 

J  J'-lhl  J-|h|-|i2l 

EE  E 


4  coshA"i^,i3,i3(/i  +  ^)  .  ^  ^  ^ 

^n,«2,»3  1  rr  7  sm  (^lOl  +  ^202  +  *3*53) 

ti=-j  t2=-./+Kj|  i3=-J+|n|+K2|  cosn/\,j,i2,i3n 


(4.12) 
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Sfi  —  {kx.n^  "I"  ky.ny  “I"  ^  —  [l  i  2,  3] 

•^1*1,12,13  ^  ^  "I"  l2^x,i2  (^l^j/.i]  *4”  ^2^1/. 12  "4*  ^3^j/,I3  )^ 


This  form  for  the  potential  function  exactly  satisfies  mass  conservation,  in  the  form 
of  the  Laplace  equation,  and  the  bottom  boundarj’  condition  for  a  locally  horizontal 
bottom.  It  allows  for  high  order  representation  of  each  of  the  stead}'  waves,  as  well 
as  the  interaction  of  each  wave  with  every  other  wave.  The  balance  of  the  solution  is 
specified  by  the  full  three  dimensional  free  surface  boundary  conditions,  Eqs.  4.4  and 
4.6. 


4.3.2  Dynamics 

The  dynamics  are  available  through  the  Bernoulli  equation  in  three  dimensions: 

^  +  +  +  +  =  0  (4.13) 

The  dynamic  pressure  is  the  difference  between  the  total  and  hydrostatic  pressure 
{pd=p-ph). 


The  Bernoulli  Constant 


In  Chapter  2,  an  explicit  and  exact  expression  for  the  Bernoulli  constant,  B,  is 
given  for  steady  two  dimensional  waves.  (Eqs.  2.14  and  2.15).  A  similar  expression 
can  be  established  for  unsteady  three  dimensional  waves. 

The  Bernoulli  equation  is  applied  at  the  bed: 


,  ^  /  2  I  2\  I 


(4.14) 


where  the  subscript,  b  indicates  the  value  at  the  bed. 

The  two  dimensional  approach  is  based  on  analysis  first  presented  by  Longuet- 
Higgins  (1975).  In  that  work,  the  analysis  was  applied  to  steady  waves.  With  steady 
waves,  the  average  over  either  time  or  space  of  the  pressure  at  the  bed  is  the  hydro¬ 
static  pressure.  With  unsteady  waves,  the  pressure  at  the  bed  averaged  over  space 
at  any  given  time,  or  averaged  over  time  at  any  given  location,  is  not  necessarily  the 
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hydrostatic  pressure.  If,  however,  an  average  is  taken  over  both  time  and  horizontal 
space,  the  result  must  be  the  hydrostatic  pressure. 


Pb  =  pg{n  -  Zb)  (4.15) 

with  the  double  overline  indicating  an  average  over  time  and  horizontal  space.  When 
the  Bernoulli  equation  at  the  bed  is  averaged  over  time  and  horizontal  space,  the 
average  time  gradient  of  the  potential  function  is  zero,  resulting  in  a  simple  expression 
for  the  Bernoulli  constant: 

B  =  ^(“6  +  ^b)  +  9^  (4.16) 

For  the  potential  function  given  by  Eq.  4.10,  and  2:  =  0  at  the  mean  water  level 
(77),  B  becomes 


B=\(Ul  +  Ul)  +  \Y, 


n=-J  «v=-J+E^=}  M 


giving  a  complete  expression  for  5  as  a  function  of  the  parameters  of  the  potential 
function. 


4*4  A  Local  Two-Intersecting- Wave  Theory 

While  a  large  number  of  intersecting  waves  could  capture  a  sea  state  of  virtually 
any  complexity,  it  would  be  difficult  to  distinguish  between  the  effects  of  each  in¬ 
dividual  wave.  It  is  important  to  remember  that  the  familiar  directional  spectrum 
description  of  a  sea  state  is  an  integral  description.  While  there  are  many  different 
frequency  and  direction  modes  represented  in  the  spectrum,  they  do  not  necessarily 
all  play  a  significant  rule  in  the  kinematics  the  entire  time.  In  fact,  when  observing 
irregular  seas,  it  is  often  the  case  that  at  any  given  time,  there  appears  to  be  a  single 
dominant  wave  of  a  particular  frequency  and  direction.  Over  time,  there  is  a  series  of 
such  dominant  waves,  each  of  a  different  frequency  and  direction.  The  integral  effect 
of  this  process  is  a  broad  directional  spectrum,  but  if  time  is  separated  into  individual 
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small  windows,  in  each  window  there  may  be  only  one  or  two  free  waves  dominating 
the  motion. 

Taking  advantage  of  the  local  nature  of  the  approach,  only  one  or  two  intersect¬ 
ing  free  waves  are  considered  in  each  window.  While  unlikel}*  to  capture  the  full 
complexity  of  a  broadly  directional  sea.  this  approach  should  be  able  to  capture  the 
dominant  modes  in  each  window.  The  direction,  frequency,  and  amplitude  of  the 
dominant  modes  will  vary  substantially  from  window  to  window  as  they  do  in  the 
two  dimensional  approach,  having  an  integral  effect  that  includes  a  large  variation  in 
directionality. 

The  two  wave  method  is  expected  to  be  particularly  effective  for  the  case  of 
standing  waves  or  short  crested  waves,  as  would  be  found  near  a  reflecting  surface 
when  the  incident  wave  field  is  almost  unidirectional,  as  it  often  is  in  shallow  water. 

Expanding  Eq.  4.11  to  fourth  order,  the  potential  function  for  two  intersecting 
waves  is: 

20 

4>{x,  y,  t)  —  UrT  -f  Uyy  -f  ^2  0C(A',7j,  A'.-z)  sin  cr,  (4.18) 

«=i 


where: 


(C(A',A,  Kiz) 


cosh  Ki{h  +  z) 
cosh  A', A 


•Si  =  -  uJit  4-  Qi) 

S2  =  {kx,23'  +  ky^2y  —  wzf  +  02) 


a; 


At  first  order: 


or,  =  (5,) 


(72  =  {S2) 


<T3  =  {2Si)  (r,  =  {2S2) 

(75  =  {Si  -f  ^2)  <7s  =  (5,  —  S2) 


At  second  order: 
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At  third  order: 


=  (35i) 

<79  =  (25’i  +  S2) 
<^11  =  +  2^2) 


(Tg  =  (352) 
o’lo  =  (25i  —  52) 
0’12  =  (Si  -  2S2) 


And  at  fourth  order: 


<?-i3  =  (45i) 

<^15  =  (35i  +  52) 
•^17  =  (25i  +  252) 
<719  =  (5i  +  352) 


<Ti4  =  (452) 

<716  =  (35l  —  52) 
<7i8  =  (25i-252) 
<720  =  (5i  —  352) 


<7i, <72, <73, 0-4, <77,(78,  <7i3,  and  <Ti4  are  self-wave  interactions.  The  balance  of  the  <7,  are 
wave-wave  interactions. 

Because  the  phases  are  arguments  of  the  sine  in  the  potential  function,  a  sim¬ 
ple  expansion  of  Eq.  4.11  results  in  redundant  terms  that  have  been  removed.  For 
example,  if  there  are  two  components: 


Bi  sin(5i  -  52)  and  B2sin(S2  -  Si)  (4.19) 

Because  sine  is  an  odd  function: 

B2  sin(52  -  5i)  =  -B2  sm(Si  -  52)  (4.20) 

the  two  components  can  be  combined  into  a  single  component: 

(Bi  -  B2)  sin(5i  -  S2)  =  Ai  sin(5i  -  52)  (4.21) 

If  both  components  were  included,  the  effects  of  Bi  and  B2  would  be  indistinguishable 
in  the  optimization;  they  must  be  combined  into  the  single  coefficient,  A,-. 

This  form  for  the  potential  function  results  in  28  unknowns  at  fourth  order  (20A,- 
plus  kj;,  ky,  a>,  and  a  for  each  of  the  two  intersecting  waves)  that  must  be  found 
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for  each  window  in  time.  At  least  28  independent  equations  must  be  defined  by 
applying  a  combination  of  each  of  the  free  surface  boundary  conditions  (Eq.  4.4  and 
4.6)  and  the  observational  equations  at  a  number  of  nodes  throughout  each  window. 
The  specific  combination  of  these  equations  will  be  determined  by  the  layout  of  the 
instruments. 

4.4.1  Kinematics 

In  order  to  apply  the  free  surface  boundary  conditions,  the  full  kinematics  must 
be  known.  While  the  kinematics  are  completely  specified  by  the  potential  function, 
Eq.  4.18,  some  of  the  algebra  may  not  be  obvious.  The  full  expressions  are  provided 
here. 


II 

H 

II 

f/x  +  Ai^  CC{h'ih,  Kiz)  cos  (T, 

«=i 

(4.22) 

v{x,y,z,t)  =  = 

/  o 

Uy  +  Yi  A~  0C( A',/i,  Kiz)  cos  a, 

1=1  ^ 

(4.23) 

te(i,y,z,t)  = 

/ 

=  ^AA',9C(A',/i,A'.2)cos<t, 

1=1 

(4.24) 

YAi^O:^iKih,Kiz)  cos  a, 

1=1 

(4.25) 

(4.26) 
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where: 


(£{Kih,  Kiz)  = 
^{Kih,Kiz)  = 


cosh  Ki{h-\-  z) 
cosh  Kih 
sinh  z) 

cosh  Kih 


Ki  = 


(Ti  =  {mSi  +  n^a) 
dcTi 

Qy  ”  “I" 


dai 

dx 

dai 

dt 


—  {pxkx^X  "I"  'nkx^2) 


=  —(mui  +  nL02) 


The  solution  details  for  a  given  set  of  measurements  are  problem  specific.  These 
will  be  provided  for  an  array  of  water  surface  measurements  in  Chapter  5  and  for  an 
array  of  pressure  measurements  in  Chapter  6. 
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Chapter  5 

Array  of  Water  Surface  Traces 


Water  surface  measurements  are  difficult  to  obtain  and  not  very  commonly  used 
for  field  measurements,  but  they  are  routinely  used  in  the  laboratory.  Surface  piercing 
wave  gauges  are  the  most  commonlj'  and  easily  deployed  of  the  methods  for  measuring 
waves  in  the  lab.  Unfortunately,  even  large  arrays  of  such  wave  gauges  still  only  give 
information  about  the  fluctuations  of  the  water  surface.  The  underlying  kinematics 
must  still  be  predicted  with  a  wave  theory.  The  following  is  a  method  for  determining 
the  kinematics  of  waves  in  the  region  of  an  array  of  water  surface  measurements.  It 
is  an  extension  of  the  LFI  method  presented  in  the  previous  chapters,  expanded  to 
determine  the  directional  structure  of  the  wave  field.  The  method  is  fully  nonlinear, 
and  results  in  a  complete  prediction  for  the  full  kinematics  in  the  vicinity  of  an  array 
of  water  surface  measurements,  throughout  the  depth  of  the  water  column. 


5.1  Formulation  of  Solution 

As  described  in  more  detail  in  Chapter  4,  the  flow  is  assumed  to  be  irrotational  and 
incompressible,  with  a  potential  function  that  represents  either  a  single  directional 
wave,  or  two  separate  intersecting  waves. 

When  the  measurement  is  taken  in  a  location  that  is  far  from  reflecting  surfaces, 
it  can  be  effective  and  straightforward  to  assume  that  the  wave  field  can  locally 
be  defined  as  a  segment  of  a  single  steady  wave.  This  is  quite  similar  to  the  LFI 
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method  for  a  point  water  surface  record  (Sobey  1992),  but  with  the  direction  of  the 
wave  determined.  The  varying  directional  trend  of  the  sea  state  is  accommodated 


by  determining  the  direction  in  each  window  separately,  allowing  for  a  fluctuation  of 
the  wave  direction  with  time,  from  window  to  window.  In  this  case,  the  potential 
function  is  Eq.  4.10  reduced  to  a  single  wave: 

<j){x,  z,  t)  =  UxX  +  Uyy  +  ^  Aj — sin  +  kyy  -Lot  +  cx)  (5.1) 


j=i 


where  Ux  and  Uy  are  the  components  of  the  known  depth  uniform  Eulerian  current 
in  the  x  and  y  directions,  h  is  the  mean  water  depth,  J  is  the  truncation  order  of  the 
Fourier  series,  Aj  are  the  Fourier  coefficients,  a;  is  the  local  fundamental  frequency, 
kx  and  ky  are  the  components  of  the  local  wave  number  in  the  x  and  y  directions, 
and  A  is  the  magnitude  of  the  local  wave  number* 

When  an  array  of  water  surface  gauges  is  placed  near  a  reflecting  surface,  such 
as  a  sea  wall,  the  resulting  sea  state  is  likely  to  contain  simultaneous  components  in 
different  directions,  such  as  in  a  standing  wave  or  short  crested  sea.  This  effect  can 
not  be  captured  with  a  potential  function  representing  a  single  progressive  wave.  It  is 
possible,  however,  to  capture  this  type  of  sea  with  a  potential  function  representing 
two  intersecting  waves.  The  potential  function  for  two  intersecting  waves  is  Eq.  4.18, 
which  is  repeated  here: 

20 

'S:,  t)  =  UxX  -f  Uyy  -f  ^  Ai  GC( A'jA,  Kiz)  sin  ai  (5.2) 

«=i 


where: 


(£{Kih,Kiz) 


cosh  Ki{h  +  z) 
cosh  Kih 


51  =  {kx,ix  +  ky^iy  —  uit  +  Oi) 

52  =  {kx^2X  +  ky^2y  —  ^2^  "f"  0:2) 


Ki  = 


2 
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At  first  order: 

=  (*51)  <^2  =  (<^*2) 

At  second  order: 

<T3  =  (25i)  fT4  =  (252) 

<^5  =  (-^1  +  *^2)  <^6  =  (•^>1  “  ■S'2) 

At  third  order: 

=  (35,) 

<79  =  (25,  +  52) 

<^11  =  (5,  +  252) 

And  at  fourth  order: 

<^13  =  (45,)  <t,4  =  (452) 

<^15  =  (35,  +  52)  (7t6  =  (35,  —  52) 

<7,7  =  (25,  +  252)  <7,s  =  (25,  —  252) 

<7,9  =  (5,  +  352)  <720  =  (5,  —  352) 

It  might  be  possible  to  include  a  larger  number  of  intersecting  free  waves  in  order 
to  capture  a  broadly  directional  sea,  but  it  would  introduce  additional  complication 
in  distinguishing  the  effects  of  the  individual  waves,  and  will  not  be  considered  here. 
Both  of  these  potential  functions  satisfy  the  field  equation  (Laplace,  Eq.  4.2)  and 
the  bottom  boundary  condition  (Eq.  4.3)  exactly.  The  balance  of  the  solution  is 
determined  by  the  free  surface  boundary  conditions:  the  modified  kinematic  free 


<78  =  (352) 

<7,0  =  (25,  —  S2) 
<7,2  =  (5,  —  252) 
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surface  boundary  condition  (/*  ), 

fK 


fiv  .  o/'  I 


dw 


+  — h  uv-p^  +  UW  . 

ox  ox  Ox 

,  du  ^dv  dw 

+  uv—  +  v^—  +  vw— 
Oy  dy  dy 

du  dv  ,  dw 
+  uw-^  +  wv—  +  w 


0 


dz  '  dz  '  ~  dz 
and  the  dynamic  free  surface  boundary  condition  (/^), 

,r)  d<f)  1,9  9  9.  — 

/  = +  9V  -  B  =  0 


at 


at 


with  the  Bernoulli  constant  (B)  defined  as: 


z  =  y 


Z  =  T) 


(5.3) 


(5.4) 


B  = 


(5.5) 


The  problem  of  determining  the  kinematics  of  irregular  waves  from  a  set  of  mea¬ 
sured  water  surface  traces  is  a  mathematically  better  posed  problem  than  interpreting 
a  subsurface  record.  The  flow  is  governed  by  the  elliptic  Laplace  equation  (Eq.  4.2), 
so  that  the  solution  is  determined  by  the  boundaries.  While  the  complete  boundaries 
of  the  solution  domain  are  not  known,  the  boundary  conditions  and  location  of  the 
boundary  are  known  at  both  the  top  and  bottom  of  the  solution  domain.  The  bottom 
boundary  condition  is  well  defined,  and  the  location  of  the  water  surface  is  measured 
at  a  few  locations  in  space  and  many  points  in  time,  allowing  the  direct  application  of 
the  free  surface  boundary  conditions.  This  is  in  contrast  to  working  with  subsurface 
records,  in  which  the  location  of  the  free  surface  must  be  determined  in  order  to  apply 
the  free  surface  boundary  conditions.  The  need  for  horizontal  boundary  conditions  is 
eliminated  by  the  assumed  periodicity  of  the  chosen  potential  function.  However  the 
fundamental  wavelength(s)  and  period(s)  must  be  found  as  part  of  the  solution. 

When  working  with  a  point  measurement,  the  fundamental  frequency  is  fairly  well 
defined  by  the  time  evolution  in  the  window  chosen,  as  long  as  the  window  is  wide 
enough.  There  is  no  direct  information  available  about  the  spatial  evolution  of  the 
signal,  however,  so  the  wave  number  is  determined  only  through  the  application  of 
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the  free  surface  boundary  conditions.  An  array  of  measurements  provides  information 
about  the  spatial  evolution  of  the  signal,  helping  to  better  define  the  fundamental 
wave  number.  This  usually  results  in  faster  and  more  robust  convergence  of  the 
optimization. 

5.2  Formulation  of  the  Optimization 

The  formulation  of  the  optimization  for  the  LFI  method  as  applied  to  the  in¬ 
terpretation  of  an  array  of  water  surface  measurements  hcis  a  great  deal  in  common 
with  the  formulation  for  a  subsurface  pressure  record  (chapter  3)  and  a  point  water 
surface  measurement  (Sobey  1992).  Less  detail  is  presented  here  than  in  the  previous 
chapter,  but  the  framework  will  be  presented,  with  an  emphasis  on  the  additional 
information  necessary  for  applying  the  method  to  an  array  of  measurements. 

Observational  Equations  The  governing  equations  presented  in  the  previous  sec¬ 
tion  represent  a  free  surface  potential  flow,  with  one  or  two  components  propagating 
in  an  arbitrary  direction.  The  observ'ational  equations  are  the  equations  in  the  system 
that  force  the  solution  to  fit  the  given  measured  record.  As  the  location  of  the  water 
surface  has  been  measured,  these  are  the  free  surface  boundary  conditions  (Eqs.  5.3 
and  5.4),  applied  at  the  horizontal  location  of  each  of  the  nodes  in  the  array.  Sufficient 
independent  equations  are  defined  by  applying  the  boundary  conditions  at  a  number 
of  times  along  the  measured  records,  within  the  window  in  time  considered  (Fig.  5.1). 
The  solution  is  the  set  of  parameters  in  the  potential  function  that  result  in  the  least 
error  in  the  FSBCs. 

In  order  to  specify  the  solution,  there  must  be  at  least  as  many  independent 
equations  as  there  are  unknown  parameters  in  the  system  of  equations.  The  free 
surface  boundary  conditions  (/,^  „  and  /^.n)  applied  at  each  of  the  N  measured 
locations  and  at  M  time  samples  in  the  window,  resulting  in  2MN  independent 
equations.  In  the  single  wave  case,  there  are  A  J  unknown  parameters  sought  in 
Eq.  5.1  {kx,  ky,  uj,  q,  and  Ai ...  Aj)  in  each  wdndow,  so  that  if  2MN  >  4  +  J,  the 
solution  is  specified.  In  the  tw'o  wav'e  case,  there  are  2^2  J  +8  unknowns  in  Eq.  5.2 
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Figure  5.1:  Schematic  of  system  of  equations  in  a  window 


(2^2:?  2uj,  2a,  and  A1...A1:  10  at  first  order,  14  at  second  order,  20  at  third 
order,  and  28  at  fourth  order).  The  solution  is  specified  when  2MN  >  2J2j  +  8- 
This  system  results  in  the  following  least  squares  optimization. 

M  AT 


minimize  0(X)  =  EE/i:',  Xn^yn’)'nm,n^tm)  +  ^n?  J/n?  (^-6) 


m=l  n=l 


where: 


X  (fcj,,  ky^U)^  (X,  Aj,  . . .  ,  Aj^ 

for  the  one  wave  case,  and: 

X  =  ky^i,uji,0!i,kx,2^  ky^2i^2i  ^2?  Ai, . . .  ,  Aj) 

for  the  two  wave  case.  Xn  and  are  the  coordinates  of  the  nth  gauge,  and  T}m,n  is 
the  measured  elevation  of  the  water  surface  at  time  tm  and  gauge  n. 

As  with  the  analysis  of  a  pressure  record,  overspecification  can  be  helpful  in  accom¬ 
modating  the  measurement  errors  in  field  records.  More  than  the  minimum  number 
of  time  samples  in  the  window  may  be  required  to  define  the  shape  of  the  water 
surface  in  each  window.  This  is  less  likely  to  be  necessary  with  two  waves  than  with 
one,  as  the  number  of  unknown  parameters  is  much  larger.  It  should  be  kept  in 
mind,  however,  as  it  would  be  a  factor  for  low  order  solutions  with  a  large  number  of 
measurement  locations. 
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5.3  Computation  Methods 

The  LFI  method  for  an  array  of  water  surface  measurements  can  be  broken  down 
into  a  similar  sequence  of  steps  as  with  the  analysis  of  a  pressure  record: 

1.  Pre-processing  of  record. 

(a)  Determine  estimate  for  level  of  noise  in  the  record. 

(b)  Determine  estimate  for  magnitude  and  direction  of  Eulerian  current. 

(c)  Define  a  set  of  continuous  records  from  each  gauge  from  the  discrete  ob¬ 
servations  by  cubic  spline  interpolation. 

(d)  Specify  spacing  of  output  locations. 

(e)  Compute  mean  zero  crossing  frequencj". 

(f)  Non-dimensionalize  record  and  all  parameters. 

2.  Primary  values  of  numerical  solution  parameters  are  chosen. 

(a)  Window  width  (tq) 

(b)  Order  of  solution  (J) 

(c)  Number  of  time  samples  of  the  water  surface  records  within  each  window 

(M) 

3.  Global  solution  is  computed  on  an  entire  wave  to  provide  first  estimates  for 
local  optimization 

4.  For  each  selected  output  location,  a  window  in  the  record  is  defined,  and  the 
an  LFI  solution  is  computed. 

(a)  Initial  guess  for  the  optimization  is  determined  from  the  global  solution. 

(b)  Full  nonlinear  optimization  for  all  unknown  components  of  the  potential 
function  is  computed. 

(c)  Results  are  checked  for  spurious  solution. 
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i.  If  no  solution,  or  a  spurious  solution,  is  found  the  solution  parameters 
are  adjusted,  and  the  optimization  repeated. 

ii.  If  a  good  solution  is  found,  progress  to  the  next  window. 

5.3.1  Pre-Processing  of  Record 

Accommodating  Measurement  Error 

Measurement  error  can  be  a  major  source  of  difficulty  with  any  high  order  data 
interpretation  method.  Local  methods  can  be  especially  sensitive,  as  each  window 
solution  relies  on  the  detail  contained  in  a  small  segment  of  the  record.  In  applying 
the  LFI  method  to  a  single  pressure  trace,  Sobey  (1992)  found  it  necessary  to  apply 
a  simple  moving  average  filter  to  field  and  laboratory  data.  In  his  work,  the  primitive 
kinematic  free  surface  boundary  condition  was  used,  requiring  an  estimate  for  the 
local  gradients  in  the  water  surface.  In  the  current  work,  a  modified  version  of  the 
boundary  condition  (Eq.  5.3)  is  used  which  does  not  require  these  gradients.  This 
makes  the  method  less  sensitive  to  noise,  so  it  was  not  necessary  to  apply  any  filtering 
for  the  results  presented  here. 

If  there  is  substantial  noise  in  the  measured  record,  the  system  of  equations  can  be 
substantially  overspecified,  allowing  the  least  squares  optimization  to  accommodate 
the  noise  in  the  record.  When  this  is  possible,  it  is  preferable  to  applying  a  smoothing 
filter  to  the  record,  as  it  does  not  impose  any  assumptions  on  the  nature  of  the  record. 
However,  if  the  error  bands  are  very  large  on  the  data,  it  may  still  be  necessary  to 
apply  filtering  to  the  raw  measurements. 

The  Mean  Water  Level 

The  mean  water  depth,  h,  must  be  specified  as  part  of  the  potential  function. 
As  a  time  series  of  the  water  surface  is  provided,  it  is  a  simple  matter  to  compute  a 
mean  of  the  measured  records.  The  mean  should  be  taken  over  a  period  much  longer 
than  a  typical  wave  period,  but  short  enough  to  accommodate  changes  in  the  local 
water  level  due  to  astronomical  and  storm  tides.  In  keeping  with  the  local  nature 
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of  the  approach,  this  is  the  local  mean  water  level,  rather  than  a  global  still  water 
level,  which  might  be  different,  and  would  be  less  appropriate.  As  with  subsurface 
mecisurements,  the  precision  of  this  value  is  not  critical,  as  h  is  only  used  to  locate 
the  origin  of  the  coordinate  system  for  the  potential  function. 

The  Eulerian  Current 

The  Eulerian  current,  Ux  and  (/y,  completes  the  description  of  the  propagation 
medium  of  the  waves,  and  an  accurate  estimate  of  its  magnitude  and  direction  is 
important.  The  measured  water  surface  traces  provide  no  information  about  the 
current  field,  so  the  information  needs  to  be  provided  from  other  sources.  See  Section 
3.3.1  for  a  detailed  discussion  of  this  important  parameter. 

Spline  of  the  Records 

As  with  the  previous  chapter,  a  set  of  continuous  records  is  computed  by  cubic 
spline  interpolation  among  the  measured  points  at  each  gauge,  allowing  complete 
flexibility  in  the  choice  of  window  widths  and  location  of  samples  in  time  of  the 
records. 

Output  Locations 

The  spacing  of  the  desired  output  locations  must  be  chosen  to  determine  the 
placement  of  the  windows.  Each  window  is  computed  independently  so  there  is  no 
restriction  on  the  spacing  of  the  output  windows.  In  addition  to  selecting  the  output 
spacing  in  time,  in  can  be  useful  when  using  an  array  of  instruments  to  determine  a 
single  central  location  within  the  array  at  which  to  define  the  solution. 

Non-dimensionalization 

In  order  to  prevent  spurious  solutions  due  to  the  comparisons  of  errors  of  different 
dimensional  quantities,  all  parameters  and  variables  are  non-dimensionalized  before 
computation  with  scales  defined  by  the  same  physically  identifiable  parameters  used 
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in  Chapter  3:  the  mass  density  of  water  (p),  acceleration  of  gravity  (g),  and  the  mean 
zero  crossing  frequency  of  the  measured  record  (Uz). 

Length  scale  =  gjuil 

Time  scale  =  l/uJz  j-j 

pg^ 

Mass  scale  = 

z 

5.3.2  Optimization  Procedure 

The  nonlinear  optimization  in  the  LFI  method  for  the  analysis  of  arrays  of  water 
surface  measurements  is  very  similar  to  that  used  for  the  analysis  of  a  subsurface 
pressure  trace.  In  the  single  wave  approach,  the  solution  is  somewhat  easier.  The 
potential  function  used  by  the  single  wave  approach  is  almost  the  same  as  that  used 
with  a  point  measurement  (both  Chapter  3  and  Sobey  (1992)),  with  the  addition 
of  a  directional  component  to  the  wave  number.  This  results  in  a  single  additional 
unknown,  but  the  array  of  measurements  provides  at  least  three  times  as  much  data, 
with  gauges  at  at  least  three  spatial  locations  to  specify  uniquely  the  directional 
structure  of  the  sea.  The  optimization  tends  to  converge  more  rapidly  and  robustly 
than  with  a  single  point  measurement. 

'When  applying  the  method  with  the  two  wave  potential  function,  there  are  many 
more  unknowns,  and  the  optimization  once  again  becomes  somewhat  tenuous.  As 
with  the  previous  chapter,  in  is  essential  to  identify  a  good  initial  estimate  to  reduce 
the  chances  that  the  optimization  will  converge  to  a  spurious  minimum. 

Global  Solution 

The  strength  of  local  methods  is  that  they  seek  a  single  solution  for  only  a  short 
segment  of  a  record  at  a  time.  The  downside  is  that  a  single  short  segment  often 
may  not  contain  sufficient  information  to  identify  the  directional  nature  of  the  wave 
field.  In  the  method  presented  in  the  previous  chapter  the  initial  estimate  could  be 
computed  from  the  data  in  the  current  window.  In  contrast,  when  working  with 
spatial  arrays,  it  is  necessary  to  examine  a  larger  segment  of  the  record,  for  example 
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an  entire  wave  from  crest  to  crest  or  trough  to  trough,  to  get  a  general  sense  of  the 
directional  trend  of  the  wave  field.  This  global  solution  can  then  be  refined  to  fit  a 
small  segment  very  closely.  This  approach  is  used  to  establish  the  initial  estimate  for 
the  optimization  procedure  in  each  window. 

For  the  initial  estimate  to  the  global  solution,  it  is  assumed  that  the  water  surface 
records  can  be  approximated  with  linear  wave  theory: 


T}  —  a  cos  (kx^n  "V  kyyn  —  a)  (5-8) 

for  the  single  wave  method,  and 


■q  =  di  cos  {kx^Xr^  +  k^^,ly„  -  iJitm Oi) 
+  02  cos  {kx.2Xr^  +  ky^^yn  -  +  <>2) 


(5.9) 


for  the  two  wave  method,  where  On  =  AnU^fg,  and  An  is  the  amplitude  of  the  linear 
potential  function  of  the  nth  wave. 


Directional  Trend  The  first  step  is  to  determine  the  directional  trend  of  the  wave 
field.  This  determination  is  accomplished  by  examining  the  gradients  of  the  water 
surface  throughout  the  segment  considered.  Estimates  for  the  spatial  gradients  and 
the  elevation  of  the  water  surface  at  the  center  of  the  array  are  computed  by  finite 
difference  approximations.  The  water  surface  is  expanded  in  a  two  dimensional  Taylor 
series  about  the  center  of  the  array: 

q{x,  y)  =  q{xc,  yc)  +  (r  -  +  (y  -  y,)^  +  . . .  (5.10) 

where  Xc  and  y^  are  the  coordinates  of  the  center  of  the  array.  This  expansion  is 
written  for  each  of  the  N  measured  gauges,  at  each  of  the  M  points  in  time,  resulting 
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where  and  are  the  coordinates  of  the  nth  gauge  in  the  array,  7/cm?  and 
are  the  water  surface  and  gradients  of  the  water  surface  at  the  center  of  the 
array  at  time,  tm,  and  rin,m  is  the  water  surface  elevation  at  gauge  n  and  time  tm-  If 
there  are  three  gauges,  the  gradients  are  uniquely  specified.  If  there  are  more,  the 
system  is  solved  in  the  least  squared  sense. 

The  water  surface  is  traveling  either  toward  or  away  from  the  direction  of  the 
water  surface  slope  depending  on  whether  it  is  moving  up  or  down  at  the  time.  The 
direction  is  thus  determined  by  the  spatial  gradient  of  the  water  surface,  and  the  sign 
of  the  time  gradient,  yielding  a  set  complex  direction  vectors: 


(5.12) 


A  cubic  spline  of  the  water  surface  at  the  center  of  the  array  (rjc)  would  provide 
an  explicit  piecewise  polynomial  expression  that  could  be  directly  differentiated  to 
obtain  the  gradients  in  time  of  the  water  surface.  These  estimates  would  be  very 
sensitive  to  measurement  errors  in  the  record.  To  obtain  more  robust  estimates  for 
the  time  gradients,  a  smoothing  cubic  smoothing  spline  (de  Boor  1978)  is  employed 
instead.  This  algorithm  provides  a  smoothing  parameter,  p,  that  can  be  set  at  any 
value  between  0  and  1,  where  p  =  0  results  in  the  linear  regression  fit,  and  p  =  1 
results  in  the  “natural”  cubic  spline.  The  smoothing  parameter  may  be  varied  to 
accommodate  varying  levels  of  noise  in  the  record.  For  the  examples  in  this  chapter, 
p  =  0.9  was  found  to  be  satisfactory. 
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For  the  single  wave  method,  it  is  expected  that  the  individual  directions  will  vary 
around  a  central  dominant  direction,  and  the  propagation  direction  estimate  is  the 
orientation  of  the  mean  of  the  complex  direction  vectors: 

e  =  angle  ^  (5.13) 

For  the  two  wave  method,  it  is  expected  that  the  propagation  directions  at  each 
point  in  time  will  vary  around  two  dominant  directions.  In  the  case  of  a  standing 
wave,  these  two  directions  are  clearly  defined.  In  a  more  complex  sea,  it  is  not  so 
straightforward.  To  separate  the  two  dominant  propagation  directions,  the  set  of 
direction  vectors,  is  divided  into  two  sets.  The  two  dominant  directions  are 
defined  as  the  angles  of  the  means  of  the  sets  of  direction  vectors  within  tt  greater 
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than  and  less  than  the  mean  direction  (see  Fig.  5.2): 

/  1  \ 

Oi  =  angle  (  ^  1  where:  0  <  angle(^^)  <  ^  +  tt 

02  =  angle  I  —  ^  where:  0  -  tt  <  a.ngle0^)  <  0 


(5.14) 


Frequency  Once  the  dominant  directions  have  been  identified,  the  approximate 
frequency  must  be  determined.  This  is  accomplished  by  examining  the  water  surface 
at  the  center  of  the  array,  as  interpolated  by  the  finite  difference  approximation 
described  above  (Eq.  5.11).  The  method  used  is  identical  to  that  used  in  the  two 
dimensional  method,  (Ch.  3): 


= 


/ 


1  d^TJe^rn 
ric,m  dt^ 


(5.15) 


where  is  computed  from  the  smoothed  spline  used  for  to  compute  for  the 

direction  estimate.  If  the  wave  field  is  unimodal  it  is  expected  that  there  will  be  a 
single  dominant  local  frequency.  The  calculated  frequencies  at  each  time  step  will  be 
similar  in  magnitude,  and  the  mean  over  the  record  is  used  as  the  first  estimate  of  the 
frequency  for  the  single  wave.  In  the  two  wave  method,  it  is  assumed  that  the  bimodal 
sea  is  the  result  of  reflection,  so  that  the  frequency  of  the  incident  and  reflected  modes 
should  be  the  same,  and  the  mean  frequency  is  used  as  the  first  estimate  for  both 
waves. 


Wave  Numbers  Once  the  frequency  is  known,  the  wave  numbers  are  estimated 
from  the  linear  dispersion  relation. 

{u  —  KjiUn)  =  tanh  A„/i,  =  A„  cos^„,  =  A„  sin  (5.16) 

where  Un  is  the  component  of  the  Eulerian  current  in  the  direction  of  the  nth  wave. 

u„  =  cos  ^tan-' 


(5.17) 
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Amplitudes  and  Phases  The  amplitudes  and  phases  of  a  particular  record  can 
be  found  by  rearranging  the  equations  as  a  linear  least  squares  problem  by  separating 
the  cosine  and  sine  components  as  was  done  in  chapter  3: 


r]c  =  a  cos  {k^x  +  k^y  —  wt  +  a) 

=  6  cos  {kj-T  +  kyy  —  ojt)  +  csin  {kj-x  +  kyy  —  u^t) 
a  = 

a  =  arctan  {—c/b) 
for  the  single  wave  method,  and 


(5.18) 


r)c  =  Oi  cos  {kj^^ix  +  ky^iy  -  wt  +  Oj)  +  02  cos  {kr,2r  +  ky,2y  -  u}t  +  02) 

=  61  cos  +  ky  iy  -  uit)  +  Ci  sin  (/tx.iX  +  ky  iy  -  oJt) 

+62  cos  (A:x,2ir  +  ky,2y  -  uJt)  +  C2  sin  {kj._2X  +  ky_2y  -  u)t) 

fli  =  \Jb\-{-c\  (5.19) 


Oi  =  arctan  (— Cl /61) 
02  =  arctan  ( — C2  762 ) 


for  the  two  wave  method.  The  system  is  determined  if  i]c{i)  is  defined  at  at  least  two 
points  in  time  for  the  single  wave  method,  and  at  four  points  in  time  for  the  two  wave 
method.  The  system  is  solved  in  the  least  squared  sense  in  the  case  of  more  points. 


Refining  the  Linear  Estimates  These  procedures  result  in  very  rough  estimates 
for  the  parameters  of  two  intersecting  linear  waves.  The  estimates  are  then  refined 
by  optimizing  for  a  best  linear  wave  theory  fit  to  the  record: 

K  M 

minimize  0(X)  =  Y.  E  "  (5-20) 

n=l  m=l 

where: 


(X.  Xfi ,  f m  )  —  O  cos  (^x^n  "1"  ®) 
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{^kxy  ky^  Of,  cCj 

for  the  single  wave  method,  and: 

/  (X,  Vn^im)  ~  0,1  cos  (kx,\^n  “h  ^y,iyn  “  "h  Oi) 

+  0,2  COS  {kx,2^n  “h  ky^2yn  —  +  O2) 


^  (^ar,l7  ky^i^  Oi,  kx,2')  ^j/,27  ^1^2,  ^17^2) 

=  y/kln+kl^ 


for  the  two  wave  method.  The  frequencies  are  determined  from  the  linear  dispersion 
relation: 


u;„  =  v/s'^ntanhfcn/i  +  (5.21) 

where  C/„  is  defined  by  Eq.  5.17  This  optimization  results  in  a  linear  estimate  for  the 
water  surface  that  fits  the  measured  records  most  closely  in  the  segment  considered. 

This  procedure  has  been  performed  on  a  segment  of  the  records  large  enough  to 
resolve  the  directional  structure  of  the  wave  field,  usually  a  complete  wave  from  crest 
to  following  crest,  or  trough  to  following  trough.  The  final  step  in  computing  the 
initial  estimates  for  the  final  window-by-window  optimization  is  to  compute  a  full 
order  global  solution  to  this  larger  segment  of  the  record.  The  initial  parameters 
for  this  full  order  global  optimization  are  provided  by  the  computed  linear  fit  to  the 
water  surface  records,  with  the  amplitudes  adjusted  for  the  potential  function: 


The  higher  order  Fourier  amplitudes  are  all  initially  set  to  zero. 

Using  these  linear  wave  theory  estimates  as  the  first  guess,  the  full  order  optimiza¬ 
tion,  (Eq.  5.6),  is  computed  to  determine  the  best  full  order  fit  to  a  global  segment  of 
the  record.  The  parameters  of  the  potential  function  computed  by  this  optimization 
are  then  used  as  the  initial  estimate  for  the  final  optimization  in  each  defined  small 
window  in  time. 
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Phase  Shift  The  phase  parameters,  oj  and  02  are  set  for  the  phase  of  the  global 
segment  of  the  record  used.  The  local  time  in  each  window  is  set  so  that  t  =  0  at 
the  center  of  the  window.  The  initial  estimate  for  the  phase  in  each  window  must  be 
shifted  to  accommodate  the  change  in  the  time  reference  frame. 

O'n  ~  -f-  On, global  (5.23) 

where  At  is  the  difference  between  the  time  in  the  two  reference  frames. 

Nonlinear  Optimization  in  Windows 

Once  the  global  solution  is  computed,  it  is  used  as  the  first  guess  for  the  parameters 
in  each  window,  and  these  parameters  are  refined  by  solving  Eq.  5.6  with  a  standard 
nonlinear  optimization  routine.  For  the  results  given,  the  Levenberg-Marquardt  algo¬ 
rithm  was  used  as  implemented  by  the  MaTLAB  Optimization  Tool  Box  (Grace  1992). 
If  the  optimization  successfully  converges  to  a  minimum,  the  solution  is  checked  to 
see  if  it  is  a  clearly  spurious  solution.  Spurious  solutions  can  be  identified  by  the  same 
criteria  used  in  chapter  3: 

•  Very  large  or  systematically  variable  errors. 

•  First  order  amplitude  smaller  than  one  of  the  higher  order  amplitudes. 

•  Unrealistically  large  or  small  frequency  or  wave  number. 

•  Large  discontinuity  between  the  windows  in  the  predicted  kinematics. 

As  in  chapter  3,  it  is  unusual  for  the  routine  to  converge  to  a  spurious  solution.  It  is 
far  more  common  for  the  routine  not  to  converge  at  all. 

If  no  solution  or  a  spurious  solution  is  found,  it  is  necessary  to  revise  the  pa¬ 
rameters  of  the  solution  as  was  discussed  in  section  3.3.2.  These  revisions  include 
increasing  the  width  of  the  window,  and  if  that  is  not  successful,  decrezising  the  order 
of  the  solution.  If  neither  of  these  adjustments  result  in  an  acceptable  solution,  the 
window  is  skipped,  and  future  analysis  must  be  interpolated  through  that  point.  As 
with  the  analysis  of  a  point  pressure  measurement,  these  adjustments  are  most  likely 
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to  be  needed  in  the  long,  flat  troughs  of  shallow  water  waves,  or  near  zero  crossings 
in  the  record.  The  additional  data  provided  by  the  array  of  measurements,  and  the 
fact  that  the  elevation  of  the  water  surface  is  measured,  provides  for  a  much  more  ro¬ 
bust  optimization  than  with  the  subsurface  pressure  measurement.  As  a  result,  these 
adjustments  are  necessary  less  often  with  an  array  of  water  surface  measurements. 

Locating  the  Water  Surface  at  the  Array  Center 

Once  the  solution  is  found,  the  potential  function,  and  thus  the  complete  kine¬ 
matics  in  the  immediate  neighborhood  of  the  array,  are  defined.  The  solution  is  likely 
to  be  most  accurate  at  the  center  of  the  array,  and  it  is  often  convenient  to  have  a 
solution  at  a  single  point,  so  the  water  surface  at  the  center  of  the  array  must  be 
found.  This  is  accomplished  by  setting  up  another  optimization  problem  with  the 
elevation  of  the  water  surface  at  a  few  nodes  in  time  throughout  the  window  as  the 
unknowns. 

The  free  surface  boundary  conditions  (Eqs.  5.4  and  5.3)  are  applied  at  the  hori¬ 
zontal  location  of  the  center  of  the  array,  at  M  points  in  time  throughout  the  solution 
window.  At  each  point  in  time,  the  only  unknown  is  the  elevation  of  the  water  sur¬ 
face,  and  the  two  boundary  conditions  provide  two  independent  equations.  The  water 
surface  is  defined  as  that  location  that  results  in  the  least  error  in  the  FSBCs,  in  the 
least  squared  sense.  Each  point  in  time  is  independent,  but  the  system  can  be  set  up 
to  solve  for  a  number  of  points  at  once.  Enough  points  should  be  found  to  specify 
the  shape  of  the  water  surface  throughout  the  window.  For  the  results  given  here, 
six  points  were  used  in  each  window  (see  Figs.  5.6  through  5.23).  The  mean  of  the 
measured  water  surface  at  the  nodes  provides  a  good  first  estimate.  This  system 
consistently  and  rapidly  converges  to  a  solution. 
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Figure  5.3:  Layout  of  water  surface  array 

5.4  Theoretical  Records 

5.4.1  Single  Wave  Records 

The  following  results  are  from  “measured”  data  generated  by  Fourier  steady  wave 
theory  (Sobey  1989).  Fourier  theory  provides  a  near-exact  solution  for  steady  waves 
that  provides  the  complete  kinematics.  This  approach  provides  a  complete  set  of  data 
to  compare  with  the  results  of  the  LFI  method,  without  the  complications  introduced 
by  the  inevitable  errors  of  data  collected  in  the  field  or  the  laboratory.  Theoretical 
records  were  used  also  because  field  or  laboratory  data  that  included  a  full  set  of 
measured  kinematics  to  compare  wdth  the  results  are  not  available. 

A  triangular  array  of  water  surface  measurements  w’as  used,  as  indicated  in  Fig. 
5.3.  The  array  is  an  equilateral  triangle  w’ith  the  same  dimensions  as  the  DWG-1 
pressure  array  (Howell  1992).  Three  measurements  were  chosen,  as  three  is  the  min¬ 
imum  number  required  to  provide  directional  information.  Additional  measurements 
would  provide  overspecification,  and  can  easily  be  accommodated  in  the  formulation. 
With  actual  field  data,  additional  measurements  are  recommended,  as  increasing  the 
number  of  instruments  would  provide  redundancy  in  the  case  of  instrument  failure, 
as  well  as  helping  to  accommodate  measurement  error. 


Figure  5.4:  LFI  predictions  ( J  =  3)  and  exact  kinematics  at  the  center  of  the 
at  the  predicted  water  surface  for  a  steady  deep  water  wave 
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Fig.  5.4  is  a  steady  deep  water  wave  generated  by  10th  order  Fourier  theory  with 
the  following  parameters:  wave  height  =  10  m,  period  =  10  seconds,  water  depth  = 
100  m,  zero  Eulerian  current,  and  direction  of  travel  10  degrees  from  the  x-axis.  This 
is  a  fairly  large  wave  in  deep  water.  The  window  width  is  1  sec,  1/10  the  zero  crossing 
period,  and  the  LFI  solution  is  computed  to  third  order  (/  =  3).  The  LFI  method  has 
captured  the  location  of  the  water  surface  and  the  kinematics  at  the  surface  essentially 
exactly.  While  linear  wave  theory  might  do  an  adequate  job  of  approximating  much 
of  the  kinematics  of  a  deep  water  steady  wave  like  this,  it  is  important  to  remember 
that  each  of  the  points  in  Fig.  5.4  was  computed  from  a  small  segment  of  the  record 
surrounding  that  point.  In  this  c«ise,  the  window  width  is  1/10  of  the  zero  crossing 
period,  or  Is.  The  local  nature  of  this  method  extends  its  applicability  to  irregular 
wave  records. 

Fig.  5.5  is  a  steady  shallow  wave  generated  by  18th  order  Fourier  theory  with  the 
following  parameters:  wave  height  =  3  m,  period  =  10  seconds,  water  depth  =  5 
m,  zero  Eulerian  current,  and  direction  of  travel  10  degrees  from  the  x-axis.  This  is 
a  fairly  extreme  wave  in  shallow  water.  The  window  width  is  1  sec.,  1/10  the  zero 
crossing  period,  and  the  LFI  solution  is  computed  to  third  order.  As  with  the  deep 
water  wave,  the  LFI  method  has  captured  the  location  of  the  water  surface  and  the 
kinematics  at  the  surface  essentially  exactly,  including  the  pronounced  sharp  crest. 

Choice  of  Order 

In  order  to  determine  the  order  necessary  to  accurately  capture  the  kinematics 
of  measured  waves,  it  is  particularly  useful  to  examine  a  window  near  the  crest  of  a 
wave.  The  crest  is  usually  the  region  that  requires  the  highest  order  solution.  This  is 
particularly  true  for  shallow  water  waves,  but  higher  order  wave  theory  in  all  depths 
of  water  indicates  that,  as  the  w'ave  height  increases,  the  crest  tends  to  get  sharper, 
and  the  trough  flatter.  Capturing  this  sharp  crest  requires  a  high  order  solution. 

Deep  Water  A  window  near  a  crest  of  the  deep  water  steady  wave  given  in  Fig  5.4 
has  been  computed  at  orders  1  through  4.  The  results  in  that  window  are  given  in 
Fig.  5.6  through  5.13.  Figures  5.6,  5.8,  5.10,  5.12  are  the  non-dimensional  errors  in 
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the  free  surface  boundary  conditions  (Eqs.  5.3  and  5.4)  at  each  of  the  measured  nodes. 
These  are  the  data  that  would  be  analyzed  in  a  practical  situation.  Figures  5.7,  5.9, 
5.11,  5.13  are  a  comparison  between  the  predicted  and  actual  values  for  the  water 
surface  and  the  velocities  at  the  surface  at  the  center  of  the  array.  The  actual  values 
were  computed  using  Fourier  wave  theory.  The  actual  values  would  not  be  available 
for  comparison  when  analyzing  field  records. 

The  first  order  computation  results  in  errors  in  the  free  surface  boundary  condi¬ 
tions  of  only  order  10"®  (Fig.  5.6)  as  well  as  very  accurately  predicting  the  velocities 
at  the  water  surface  at  the  center  of  the  array  (Fig.  5.7).  It  is  a  surprisingly  accu¬ 
rate  first  order  solution.  This  is  because  of  the  local  nature  of  the  method.  When 
a  single  first  order  solution  is  used  to  capture  the  entire  wave,  the  error  is  larger,  of 
order  10"^.  It  also  should  be  noted  that  this  first  order  solution  is  not  the  same  as  a 
linear  wave  theory  solution,  even  locally.  The  full  nonlinear  boundary  conditions  are 
preserved,  and  the  frequency  and  wave  number  are  free  to  vary,  and  are  not  bound 
by  a  dispersion  relationship.  For  steady  waves  in  deep  water,  linear  wave  theory  is 
fairly  accurate.  Linear  theory  is  not,  however,  a  local  solution,  and  is  not  directly 
applicable  to  irregular  waves. 

At  second  order,  the  free  surface  boundary  condition  errors  arc  smaller,  of  order 
10"®,  and  the  velocities  at  the  surface  match  the  Fourier  solution  visually  perfectly. 
At  third  and  fourth  order,  the  errors  in  the  free  surface  boundary  conditions  continue 
to  decline,  and  the  water  surface  velocities  continue  to  match  the  Fourier  solution 
well.  In  deep  water,  for  waves  of  this  height,  second  order  is  more  than  adequate  to 
capture  the  surface  kinematics  of  this  wave.  Higher  order  solutions  are  likely  to  be 
necessary  to  capture  the  irregularity  of  field  records,  even  in  deep  water. 

Choice  of  order  is  dictated  by  the  desired  accuracy  of  the  solution,  and  by  the  ease 
of  convergence  to  the  solution.  As  there  are  more  free  parameters  at  higher  order, 
a  higher  order  solution  will  always  have  smaller  errors  in  the  free  surface  boundary 
conditions.  In  the  case  of  this  example,  the  solution  converged  at  all  orders  very 
quickly,  so  there  is  little  penalty  in  using  third  or  fourth  order.  With  an  irregular 
record,  in  contrast,  convergence  can  be  more  difficult,  and  it  is  occasionally  necessary 
to  resort  to  lower  order. 
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As  the  order  is  increased,  there  are  more  free  parameters,  and  care  must  be  taken 
to  include  sufficient  points  in  each  window.  As  discussed  above,  for  fairly  low  orders 
and  an  array  of  measurements,  the  limiting  factor  is  the  minimum  number  of  points 
needed  to  define  the  curvature  in  the  window,  rather  than  to  specify  the  system  of 
equations.  The  above  examples  were  computed  from  an  array  of  three  points,  and 
so  required  a  minimum  of  only  one  sample  to  specify  the  system  at  first  and  second 
order,  and  two  points  for  up  to  order  8.  When  attempted,  the  solution  would  not 
converge  with  only  one  sample.  With  two  samples,  the  optimization  algorithm  found 
a  reasonable  solution,  but  this  small  number  of  samples  would  not  define  the  shape  of 
the  water  surface  adequately  if  it  were  part  of  an  irregular  record,  and  so  three  points 
were  used  in  each  window  for  all  orders.  In  general,  a  minimum  of  three  points  should 
be  used,  and  more  may  be  necessary  to  accommodate  a  highly  irregular  profile,  or 
the  inevitable  measurement  error  in  a  field  record. 

In  the  case  of  theoretically  generated  records  the  need  for  more  sampling  of  points 
poses  no  problem,  but  with  field  records,  there  are  limitations  as  to  the  spacing  of 
the  sampled  points.  In  order  to  free  up  the  spacing  of  points  for  the  LFI  method, 
points  are  sampled  from  a  cubic  spline  interpolation  of  the  actual  record.  This  allows 
the  points  to  be  sampled  anywhere  within  the  window.  While  computationally  it  is 
possible  to  sample  as  many  points  as  necessary  in  a  small  window,  if  that  window  is, 
in  fact,  defined  by  only  a  couple  of  actual  data  points,  it  is  not  appropriate  to  try 
to  fit  a  high  order  solution  to  a  segment  defined  by  only  a  few  observational  points. 
In  order  to  include  sufficient  actual  data  points  to  justify  the  increased  order,  the 
window  must  be  increased  in  size.  While  increasing  the  size  of  the  window  permits  a 
higher  order  solution,  it  also  compromises  the  local  nature  of  the  method.  The  goal 
of  the  LFI  method  is  for  the  solution  to  be  as  local  as  possible,  which  is  achieved  by 
selecting  as  small  a  window  as  possible  at  fairly  low  order. 

Shallow  Water  A  window  near  a  crest  of  the  shallow  water  steady  wave  given 
in  5.5  has  been  computed  at  orders  1  through  5.  The  results  are  given  in  Fig.  5.14 
through  5.23.  These  figures  are  analogous  to  those  previously  discussed,  but  on  a 
shallow  water  wave. 


-0.4  -0.3  -0.2  -0.1  0  0.1  0.2  0.3  0.4 

Figure  5.14;  Free  surface  boundary  condition  errors  for  a  shallow 
water  w'ave  at  order  1 


Crest  of  a  shallow  water  w'ave 
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Figure  5.15;  Free  surface  and  velocities  at  the  center  of  the  array 
for  a  shallow  water  w'ave  at  order  1 


Crest  of  a  shallow  water  wave 


Figure  5.20:  Free  surface  boundary  condition  errors  for  a  shallow 
water  wave  at  order  4 


Crest  of  a  shallow  water  wave 
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Figure  5.21:  Free  surface  and  velocities  at  the  center  of  the  array 
for  a  shallow  water  wave  at  order  4 
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In  this  case,  the  first  order  computation  results  in  substantial  errors  in  the  free 
surface  boundary  conditions  of  order  10"^  (Fig.  5.14)  and  does  not  predict  the  ve¬ 
locities  at  the  water  surface  at  the  center  of  the  array  as  well  as  it  did  in  deep  water 
(Fig.  5.15),  noticeably  under-predicting  the  horizontal  velocities.  The  local  nature 
of  the  method  allows  a  first  order  solution  to  get  fairly  close,  but  steady  waves  in 
shallow  water  are  highly  nonlinear,  and  a  first  order  solution  simply  can  not  capture 
the  sharp  crest. 

At  second  order,  the  free  surface  boundary  condition  errors  are  slightly  smaller, 
but  still  of  order  lO"^,  and  still  under  predicting  the  horizontal  velocities  at  the 
surface.  At  third  and  fourth  order,  the  errors  in  the  free  surface  boundary  conditions 
continue  to  decline,  but  the  water  surface  velocities  do  not  match  the  Fourier  solution 
visually  perfectly  until  fifth  order  (Fig.  5.23).  In  shallow  water,  it  may  be  necessary 
to  use  up  to  fifth  order  to  capture  accurately  the  surface  kinematics  of  a  fairly  large 
wave.  In  the  case  of  this  example,  the  solution  converged  at  all  orders  very  quickly, 
and  there  is  little  penalty  in  using  fourth  or  fifth  order.  With  an  irregular  record, 
convergence  will  be  more  difficult,  and  it  may  sometimes  be  necessary  to  resort  to 
lower  order. 

As  previously  discussed,  the  use  of  higher  order  solutions  requires  that  more  points 
be  sampled  in  each  window.  This  is  not  likely  to  be  a  limitation  with  the  single  wave 
form  with  an  array  of  measurements,  as  it  is  likely  to  be  necessary  to  overspecify  the 
system  in  order  to  define  the  curvature  within  the  window. 

5.4.2  Records  of  Two  Intersecting  Waves 

In  order  to  develop  the  method  using  two  intersecting  waves,  theoretical  water 
surface  records  were  generated  by  Ohyama’s  fourth  order  intersecting  wave  theory 
(Ohyama,  Jeng,  and  Hsu  1995a).  Ohyama’s  method  is  a  Stokes-style  solution  for 
irrotational  intersecting  waves  that  is  accurate  to  fourth  order,  and  applicable  in  deep 
water  (see  Appendix  A).  As  the  resulting  water  surface  from  intersecting  waves  can 
be  quite  complex,  a  desired  wave  height  can  not  be  specified.  Rather,  the  amplitude 
of  the  first  order  component  of  each  of  the  intersecting  free  waves  is  specified.  The 
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first  order  components  are  the  largest  components  of  the  resulting  wave  field,  and  thus 
give  an  approximation  of  the  size  of  the  resulting  waves.  This  wave  theory  assumes 
a  zero  Eulerian  current. 

Standing  Wave  Figure  5.24  shows  the  results  of  the  method  for  a  standing  wave 
generated  by  two  identical  intersecting  waves  traveling  in  opposing  directions.  The 
parameters  of  the  waves  are:  Period  =  10  sec.,  first  order  amplitudes  =  3  m,  prop¬ 
agation  directions,  15  and  195  degrees  from  the  x  axis,  in  deep  water.  The  LFI 
method  to  third  order  (J  =  3)  finds  the  kinematics  at  the  water  surface  almost  ex¬ 
actly.  While  these  results  show  the  complete  wave,  it  is  important  to  keep  in  mind 
that  as  with  all  the  previous  results,  each  of  the  indicated  points  is  in  the  center  of 
a  separate  window,  and  was  computed  independently  of  the  other  windows.  In  this 
case,  the  standard  window  width  was  1  sec.,  or  one  tenth  the  zero-crossing  period  of 
the  wave,  with  a  third  order  potential  function  {J  =  3),  and  four  water  surface  nodes 
distributed  equidistantly  in  time  in  each  window  (A/  =  4),  at  each  gauge,  for  a  slight 
overspecification,  with  24  equations  in  20  unknowns 

Short  Crested  Wave  Figure  5.25  is  the  water  surface  at  a  point  in  time  of  the 
short  crested  wave  that  would  result  in  the  reflection  of  a  steady  wave  from  a  sea 
wall,  as  indicated  by  figure  5.26.  The  parameters  of  the  wave  are:  Period  =  10  sec., 
first  order  amplitudes  =  3  m,  propagation  directions,  30  and  150  degrees  from  the  x 
axis,  in  deep  water.  Figure  5.27  shows  the  results  of  the  method  for  this  short  crested 
wave.  The  LFI  method  to  third  order  again  finds  the  kinematics  for  this  w'av'e  almost 
exactly.  In  this  case,  the  standard  window  width  was  also  1  sec.,  or  one  tenth  the 
zero-crossing  period  of  the  wave,  with  a  third  order  potential  function  (J  =  3)  and 
four  water  surface  nodes  (M  =  4)  in  each  window,  for  a  slight  overspecification,  with 
24  equations  in  20  unknowns. 

Choice  of  Order 

As  with  the  previous  examples  with  steady  wav'es,  it  is  useful  to  examine  a  window 
near  the  crest  in  order  to  determine  the  order  necessary  to  accurately  capture  the 


Figure  5.24:  LFI  predictions  ( J  =  3)  and  exact  kinematics  at  the  center  of  the 
at  the  predicted  water  surface  for  a  standing  wave 
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Figure  5.25:  Water  surface  of  short  crested  wave  at  t=0 


y 


Figure  5.26:  Schematic  of  reflection  that  would  generate  the  short  crested  wave  given 
in  figure  5.25 


kinematics  of  a  short  crested  wave. 

A  window  near  a  crest  of  the  short  crested  wave  given  in  Fig.  5.27  has  been  com¬ 
puted  at  orders  1  through  4.  The  results  in  that  window  are  given  in  Fig.  5.28  through 
5.35.  Figures  5.28,  5.30,  5.32,  5.34  show  the  errors  in  the  free  surface  boundary  con¬ 
ditions  at  each  of  the  measured  nodes.  These  are  the  data  that  would  be  analyzed  in 
a  practical  situation.  Figs.  5.29,  5.31,  5.33,  5.35  give  a  comparison  between  the  pre¬ 
dicted  and  actual  values  for  the  water  surface  and  the  velocities  at  the  surface  at  the 
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Figure  5.28:  Free  surface  boundary  condition  errors  for  a  short 
crested  wave  at  order  1 


Crest  of  a  short  crested  wave 


Figure  5.29:  Free  surface  and  velocities  at  the  center  of  the  array 
for  a  short  crested  wave  at  order  1 


Figure  5.31:  Free  surface  and  velocities  at  the  center  of  the  array 
for  a  short  crested  wave  at  order  2 
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Figure  5.32:  Free  surface  boundary  condition  errors  for  a  short 
crested  wave  at  order  3 
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Figure  5.33:  Free  surface  and  velocities  at  the  center  of  the  array 
for  a  short  crested  wav'e  at  order  3 
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center  of  the  array.  The  theoretical  values  used  Ohyama’s  fourth  order  intersecting 
wave  theory.  The  actual  values  would  not  be  a\'ailable  for  comparison  when  analyzing 
field  records. 

It  is  clear  that  the  first  order  computation  results  in  substantial  errors  in  the 
free  surface  boundary  conditions  of  order  10"^  (Fig.  5.28),  as  well  as  significantly 
underestimating  the  horizontal  velocities  (Fig.  5.29).  As  with  the  steady  deep  water 
wave  previously  discussed,  the  first  order  solution  is  quite  reasonable.  This  is  because 
of  the  local  nature  of  the  method,  and  the  fact  that  it  is  not  a  linear  solution.  The  full 
nonlinear  boundary  conditions  are  imposed,  and  the  frequencies  and  wave  numbers 
are  free  to  vary,  and  are  not  bound  by  a  dispersion  relationship. 

At  second  order,  the  free  surface  boundary  condition  errors  are  better,  of  order 
10"^,  and  the  velocities  at  the  surface  match  the  analytical  solution  visually  perfectly. 
At  third  and  fourth  order,  the  errors  in  the  free  surface  boundary  conditions  continue 
to  decline,  and  the  water  surface  velocities  continue  to  match  the  theoretical  solution 
well. 

As  with  the  single  wave  method,  choice  of  order  is  dictated  by  the  desired  ac¬ 
curacy  of  the  solution,  and  by  the  ease  of  convergence  of  the  optimization.  For  the 
above  examples,  four  water  surface  nodes  (M  =  4)  distributed  equidistantly  in  time 
in  each  window,  at  each  gauge,  at  orders  1  through  3.  These  values  providing  an  over¬ 
specification  and  sufficient  points  to  define  the  shape  of  the  water  surface  throughout 
the  window.  Five  points  (M  =  5)  were  used  at  fourth  order,  as  four  points  provide 
only  24  equations,  and  there  are  28  parameters  to  be  found.  Five  points  provides  30 
equations  for  a  slight  over  specification.  At  fifth  order,  there  are  38  unknowns,  and 
seven  points  in  time  (M  =  7)  would  hav'e  to  be  used  in  each  window.  Sampling  this 
many  points  in  a  single  window  w'ould  require  ver\'  closely  spaced  data  or  a  wider 
window.  Another  solution  w'ould  be  to  use  an  array  w'ith  additional  measurement 
locations.  An  array  of  four  gauges,  for  example,  would  provide  eight  equations  per 
point  in  time,  and  would  allow  a  fifth  order  solution  to  be  computed  from  five  points 
(M  =  5)  per  window'.  It  would  be  advantageous  to  use  as  many  gauges  as  possible  in 
shallow  water,  where  higher  order  solutions  are  necessary. 

Unfortunately,  the  theoretical  solution  used  for  this  analysis,  while  the  best  avail- 
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able  for  intersecting  waves,  is  only  accurate  to  fourth  order,  and  only  being  applicable 
in  deep  or  transitional  water.  This  precluded  an  analysis  of  the  behavior  of  the  so¬ 
lution  in  shallow  water.  When  the  method  is  applied  in  shallow  water,  the  lessons 
learned  in  the  two  dimensional  case  should  be  applied,  and  a  higher  order  solution, 
perhaps  to  fifth  order,  should  be  used. 

It  is  also  the  case  that  the  Ohyama  solution  is  a  single  solution  that  is  accurate 
to  fourth  order  globally.  The  local  nature  of  the  LFI  method  allows  each  separate 
window  in  time  to  be  represented  by  a  unique  solution,  representing  the  entire  record 
with  far  more  accuracy  than  a  global  solution  of  the  same  order. 

Choice  of  Window  Width 

As  mentioned  in  the  previous  section,  it  is  important  to  keep  the  window  width  at 
a  minimum.  For  the  short  crested  wave  shown  in  Figure  5.27,  the  standard  window 
width  was  one  tenth  the  zero-crossing  period  (O.IT^).  It  was  necessary  to  widen  the 
window^  to  0.157^  at  three  locations  in  order  to  find  a  solution.  This  indicates  that 
the  standard  window  width  is  set  close  to  the  smallest  size  that  would  be  effective. 
In  general,  it  is  necessary  for  the  window  to  contain  enough  of  the  record  to  define 
the  local  curvature,  and  include  sufficient  data  points  to  reasonably  expect  to  define 
a  high  order  solution.  Most  often,  the  limiting  factor  will  be  the  sample  rate  of  the 
data.  For  example,  the  above  examples  are  computed  for  a  short  crested  wave  with 
a  period  of  ten  seconds.  In  those  examples,  the  window  width  of  1/10  the  period 
(Is)  was  adequate  to  capture  most  of  the  wave.  Unfortunately,  field  records  are  often 
sampled  at  only  IHz.  At  this  sampling  frequency,  and  a  one  second  window,  there 
would  be  a  maximum  of  only  two  data  points  in  each  window,  and  it  may  be  necessary 
to  widen  the  window  to  include  more  data. 

Effects  of  Array  Size 

The  examples  above  are  all  computed  from  data  sampled  by  a  fairly  small  array 
(Fig.  5.3).  The  1.6m  size  of  the  array  is  about  0.01  of  the  wavelength  of  10s  wave  in 
deep  water.  That  particular  size  was  chosen  because  it  has  the  same  dimensions  as 
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the  DWG-1  pressure  array  that  is  used  for  the  analysis  of  pressure  records  in  Ch.  6. 
A  small  array  has  the  advantage  of  having  a  small  distance  to  interpolate  to  the 
center  of  the  array,  and  is  local  in  the  sense  that  each  of  the  sensors  is  in  essentially 
the  same  part  of  the  wave.  The  disadvantage  of  a  small  array  is  that  the  differences 
between  the  measured  values  at  each  of  the  nodes  are  small,  and  thus  more  sensitive 
to  measurement  error.  With  the  theoretical  records  used  for  the  previous  analysis, 
the  method  was  not  sensitive  to  array  size,  giving  equally  accurate  results  with  arrays 
up  to  ten  times  the  size  of  the  DWG-1.  The  spacing  between  the  gauges  of  an  array 
will  determine  the  wavelengths  to  which  the  gauge  is  most  sensitive.  Window  widths 
of  approximately  one  tenth  a  typical  period  have  been  found  to  be  effective.  This 
indicates  that  a  gauge  spacing  of  about  one  tenth  of  the  expected  dominant  wave 
length  might  be  an  optimum  spacing. 


5.5  Laboratory  Measurements 

The  results  from  the  analysis  of  analytically  derived  records  given  in  the  previous 
sections  demonstrate  the  method’s  capabilities  in  a  variety  of  conditions.  They  have 
also  been  very  useful  in  helping  to  determine  the  range  of  solution  parameters  that 
must  be  adopted,  including  order  of  solution,  window  width,  and  the  number  of  time 
samples  taken  in  each  window.  The  question  still  remains,  however,  as  to  how  well 
the  method  works  with  actual  irregular  waves. 

5.5.1  Laboratory  Experiment 

The  following  results  are  taken  from  a  laboratory  experiment  performed  together 
with  Dr.  Steven  A.  Hughes  at  the  Coastal  Engineering  Research  Center,  Army  Corps 
of  Engineers  Waterways  Experiment  Station  in  Vicksburg,  MS,  in  June  of  1995.  The 
waves  were  generated  in  a  0.46m  wide  flume  by  a  programmable,  hydraulically  driven, 
piston  type  wave  board  (see  Fig.  5.36).  At  the  other  end  of  the  flume  was  a  beach 
with  1:30  slope,  and  a  single  layer  of  wav^e- absorbing  horsehair  matting.  This  beach 
has  been  shown  to  have  bulk  reflection  coefficients  that  vary  from  about  5%  for 
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Figure  5.36:  Schematic  of  wave  flume 

Hmo/h  =  0.1  up  to  about  15%  for  Hmofh  =  0.4  (where  Hmo  is  the  incident  spectral 
significant  wave  height  and  h  is  the  water  depth)  (Sultan  1992;  Sultan  and  Hughes 
1993). 

The  water  surface  was  measured  by  capacitance  wave  rods,  calibrated  with  a  cubic 
calibration  function.  The  velocity  data  were  collected  using  a  Dantek  laser  Doppler 
velocimeter  (LDV)  system  operated  in  the  backscatter  mode.  The  LDV  system  fea¬ 
tured  a  2-watt  argon-ion  laser  equipped  with  a  fiber-optic  probe  that  measures  two 
orthogonal  water  velocity  components  (horizontal  and  vertical).  Velocity  data  were 
converted  in  real  time  to  engineering  units  (m/s)  and  written  to  a  computer  file 
simultaneously  with  the  wave  rod  data. 

The  wave  rods  and  LDV  were  placed  near  the  middle  of  the  flume,  with  the  wave 
rods  arranged  in  an  equilateral  triangle  with  leg  length  of  0.17m  (see  Fig.  5.37).  The 
LDV  was  situated  to  measure  the  horizontal  and  vertical  velocities  at  the  center  of  the 
array,  at  a  variety  of  vertical  elevations.  The  flume  was  allowed  to  reach  quiescence 
in  between  runs,  and  the  waves  were  measured  for  only  a  short  time  after  starting  the 
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Figure  5.37:  Layout  of  array  in  flume 


wave  maker.  At  a  depth  of  0.35m,  the  approximate  wave  celerity  [y/^)  is  1.8m/s.  As 
the  measuring  location  is  located  9m  from  the  beach,  it  takes  about  10s  for  the  waves 
to  travel  from  the  measurement  location  to  the  beach,  and  for  the  reflected  waves 
to  travel  back  to  the  study  location.  Thus,  the  first  10s  of  the  record  are  assured  of 
being  uncontaminated  by  reflection  from  the  beach.  Using  only  the  beginning  of  the 
record  also  prevents  a  Stokes  drift  induced  return  current  to  be  established,  and  so 
the  Eulerian  current  is  taken  to  be  zero. 


5.5.2  Laboratory  Results 

Figure  5.38  shows  a  sample  wave  taken  from  an  irregular  wave  record  measured 
in  the  flume.  The  waves  were  generated  with  a  target  TMA  spectrum  (Bouws  et  al. 
1985)  with  a  peak  period  of  2.5s  and  mean  wave  height  of  approximately  12cm.  The 
water  depth  was  0.35m,  with  the  velocities  measured  by  the  LDV  at  an  elevation  of 
0.05m  below  the  mean  water  level.  The  solid  lines  on  figure  5.38  are  the  horizontal  and 
vertical  velocities  as  measured  by  the  LDV.  The  circles  are  the  velocities  predicted  by 
the  LFI  method  at  the  center  of  each  computed  window.  The  solid  line  in  the  water 
surface  plot  is  the  measured  elev'ation  at  wave  rod  2,  at  the  same  x  coordinate  as  the 
center  of  the  array.  The  circles  are  the  water  surface  ele\’ations  at  the  center  of  the 
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Figure  5.39:  Free  surface  boundary  condition  errors  at  the  crest 
of  a  sample  wave  in  a  flume 


Figure  5.40:  Computed  and  measured  velocities  near  the  crest 
of  a  sample  wave  in  a  flume 
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array  as  computed  by  the  LFI  method. 

This  solution  was  computed  to  third  order,  with  a  window  width  of  1/10  the  mean 
zero  crossing  period  of  the  record  (approximately  0.2s).  Four  time  samples  were  taken 
in  each  window,  providing  substantial  overspecification  of  the  system,  and  allowing 
the  shape  of  the  water  surface  to  be  well  defined  in  each  window.  The  first  trough  is 
fairly  flat,  and  the  optimization  did  not  converge  to  a  solution  quickly  with  the  small 
window.  After  the  window  was  widened,  the  optimization  converged  quickly  to  the 
given  results. 

The  LFI  method  has  captured  the  detail  of  the  kinematics  of  this  irregular  wave 
very  well.  The  comparisons  to  the  measured  data  are  given  just  below  the  troughs  of 
the  wave.  Note  that  the  LFI  method  accurately  captured  the  secondary  hump  in  the 
vertical  velocity  (w)  near  the  2.5s  point.  The  data  used  to  compute  the  kinematics 
were  the  measured  water  surface.  The  computed  water  surface  was  found  from  the 
predicted  kinematics,  and  it  matches  the  measured  water  surface  almost  exactly.  The 
accurate  computation  of  the  water  surface  indicates  that  the  predicted  kinematics 
near  the  surface  may,  in  fact,  be  more  accurate  than  the  predictions  at  the  elevation 
where  the  kinematics  were  measured.  This  is  to  be  expected,  as  the  data  used  to 
compute  the  solution  was  measured  at  the  surface,  and  thus  the  results  should  be 
most  accurate  there.  Achieving  the  greatest  accuracy  near  the  surface  can  be  useful, 
as  the  surface  is  the  location  where  the  water  velocities  and  accelerations  are  greatest. 

Figures  5.39  and  5.40  show  the  results  of  the  computation  in  a  window  near  the 
crest  of  the  wave  given  in  Fig.  5.38.  The  errors  in  the  DFSBC  are  of  order  10"^,  and 
substantially  smaller  in  the  MKFSBC.  The  computed  and  measured  velocities  do  not 
match  exactly,  but  the  magnitude  and  trend  are  very  close.  The  results  are  similar 
in  the  other  windows. 


5.6  Discussion 

The  LFI  method  is  an  effective  and  efficient  way  to  re-create  the  kinematics  of 
irregular  waves  from  the  measurements  of  an  array  of  wave  rods.  It  is  able  to  capture 
the  kinematics  from  primarily  uni-directional  seas  as  well  as  the  seas  near  a  reflect- 
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ing  surface.  The  local  nature  of  the  approach  allows  a  nonlinear  solution  without 
prohibitive  computational  costs,  using  a  fairly  simple  form  for  the  potential  function, 
and  allowing  the  parameters  of  the  potential  function  to  change  with  time. 

The  examples  in  this  Chapter  were  all  computed  using  MaTLAB  running  under 
the  Linux  operating  system  on  an  Intel  90MHz  Pentium  processor  based  PC.  The 
shallow  water  steady  wave  (Fig  5.5)  took  the  longest  to  compute,  about  30min,  and 
85x10®  floating  point  operations  (flops).  The  steady  deep  water  wave  (Fig.  5.4)  took 
the  least  time,  lOmin.  and  15x10®  flops,  the  standing  wave  (Fig.  5.24)  took  20min. 
and  244x10®  flops,  the  short  crested  wave  (Fig.  5.27)  took  llmin.  and  93x10®  flops, 
and  the  laboratory  record  (Fig.  5.38)  took  21min.  and  83x10®  flops. 

The  examples  given  in  this  chapter  provide  guidance  as  to  the  parameters  of  a 
solution  to  be  applied  to  field  records.  Far  from  reflecting  surfaces,  using  a  potential 
function  representing  a  single  wave  is  effective.  Near  a  reflecting  surface,  a  potential 
function  representing  two  nonlinear  intersecting  waves  is  capable  of  capturing  the 
standing  or  short  crested  waves  that  are  likely  to  dev'elop. 

Window  widths  of  1/10  of  the  zero  crossing  period  are  small  enough  to  maintain 
the  local  nature  of  the  solution,  and  capture  the  detail  of  the  record,  while  being 
large  enough  to  include  the  local  trend  of  the  wave  field.  On  certain  segments  of  the 
record,  the  window  width  must  be  increased  in  order  to  include  sufficient  curvature 
in  the  record  to  find  a  solution.  The  widening  of  the  window  is  most  often  needed 
in  long,  flat  troughs  in  shallow  w’ater,  or  near  zero  crossings  of  the  record.  In  either 
case,  the  window  rarely  needs  to  be  larger  than  1/5  of  the  zero  crossing  period. 

Low  order  solutions  are  quite  adequate  for  deep  water  waves.  For  most  waves, 
second  order  is  adequate.  For  the  very  steepest  waves,  slightly  higher  order  may  be 
appropriate,  and  there  is  little  computational  penalty  in  including  the  higher  order 
terms.  In  shallow  water,  higher  order  solutions  are  necessary.  Third  order  is  adequate 
in  many  cases,  but  including  up  to  fifth  order  is  recommended  for  extreme  waves. 
When  attempting  a  high  order  solution  such  as  this  with  the  two  wave  method,  it 
may  be  necessary  to  use  arrays  with  more  than  three  gauges  in  order  to  provide 
sufficient  equations  to  specify  the  solution. 
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Chapter  6 

Array  of  Subsurface  Pressure 
Measurements 


Arrays  of  subsurface  pressure  gauges  are  installed  in  the  sea  in  order  to  capture 
directional  wave  field  data.  These  instruments  have  an  advantage  over  arrays  of 
water  surface  wave  rods,  as  they  are  mounted  below  the  surface,  where  they  are  less 
susceptible  to  damage,  by  either  the  strong  wave  forces  experienced  during  storms, 
as  well  as  vandalism  or  accidents  with  vessels. 

Unfortunately,  the  reason  they  are  less  susceptible  to  damage  from  wave  forces  is 
because  they  are  placed  under  the  surface,  often  at  the  sea  bed,  where  the  action  of 
waves  is  the  smallest.  This  leads  to  the  mathematically  ill-posed  problem  discussed 
in  chapter  3.  This  difficulty  is  somewhat  mitigated  by  the  added  information  made 
available  by  the  multiple  gauges  in  the  array,  but  still  limits  the  ability  to  determine 
the  detail  of  the  kinematics,  particularly  when  there  is  a  lot  of  high  frequency  energy 
in  the  sea  state. 


6.1  Formulation  of  Solution 

The  lessons  learned  in  the  previous  chapters  provide  for  a  fairly  straightforward 
determination  of  the  formulation  for  adapting  the  LFI  method  to  the  analysis  of 
arrays  of  pressure  gauges.  The  formulation  for  a  single  subsurface  pressure  gauge  was 
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presented  in  chapter  3,  and  for  arrays  of  water  surface  measurements  in  chapter  5. 
The  formulation  for  arrays  of  subsurface  pressure  gauges  is  very  similar  to  that  for 
arrays  of  water  surface  measurements,  with  the  addition  of  the  need  to  determine  the 
location  of  the  water  surface,  as  was  done  with  the  analysis  of  a  subsurface  pressure 
trace.  The  following  description  includes  the  required  information;  greater  detail  has 
been  given  in  previous  chapters. 

The  flow  is  assumed  to  be  irrotational  and  incompressible,  with  a  potential  func¬ 
tion  that  represents  either  a  single  directional  wave,  or  two  separate  intersecting 
waves.  The  potential  function  locally  representing  a  single  wave  is  Eq.  4.10  reduced 
to  a  single  wave: 

<f>{x,z,t)  =  UiT  +  Uyy  +  smjjkjrT  +  k^y  -  ut  +  o)  (6.1) 

A'  =  +  kl 

where  Ux  a,nd  Uy  are  the  components  of  the  known  depth  uniform  Eulerian  current 
in  the  x  and  y  directions,  h  is  the  mean  water  depth,  J  is  the  truncation  order  of  the 
Fourier  series,  Aj  are  the  Fourier  coefficients,  u)  is  the  local  fundamental  frequency, 
kx  and  ky  are  the  components  of  the  local  fundamental  wave  number  in  the  x  and  y 
directions,  and  K  is  the  magnitude  of  the  local  wave  number. 

The  potential  function  for  two  intersecting  waves  to  fourth  order  is  Eq.  4.18, 
repeated  here: 

20 

<f>{x,  y,  z,  t)  =  UxT  A  Uyy  -f-  ^  Ai  0C(  A'.h,  h'iz)  sin  ai  (6.2) 

t=i 


where: 
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At  first  order: 


II 

(^2  =  (52) 

At  second  order: 

=  (25i) 

<74  =  (2^2) 

<T5  =  (^i  +  5*2) 

1 

II 

At  third  order: 

<T7  =  (SSi) 

(Js  —  (3*92) 

^9  =  (25'i  +  5^2) 

1 

II 

0 

b 

<^11  =  {Si  +  2S2) 

<^12  —  ('S'l  —  2*92) 

And  at  fourth  order: 

^13  =  i^Si) 

<7i4  =  (4S2) 

<^15  =  i^Si  +  S2) 

<7i6  =  (351  ~  ^2) 

^17  =  (2*Si  +  2*92) 

<7i8  =  (25i-252) 

<719  =  (^i  +  35'2) 

<720  =  (5i  -  352) 

Both  of  these  potential  functions  solve  the  field  equation  (Laplace,  Eq.  4.2)  and  the 
bottom  boundary  condition  (Eq.  4.3)  exactly.  The  rest  of  the  solution  is  determined 
by  the  free  surface  boundary  conditions:  the  modified  free  surface  boundary  condition 


,  I  o/ 


ydu 

dx 

du 


dv 


dw 


+  uv—  +  uw- 
ox  ax 


+  uv—  + 
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+  vw 


dw 
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dy'  dy 

.  du  dv  ^dw 
+  uw—  +  wv—  +  w^—  =  0 
dz  dz  dz 
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(6.3) 


Z  =  7] 
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the  dynamic  free  surface  boundary  condition  (/®), 

/^  =  •^  +  u'^)  +  gt]  —  B  =  0  at  z  =  g  (6.4) 

the  Bernoulli  equation  (/^), 

=  ^  +  +  u^)  +  ^  B  =  0  (6.5) 

with  the  Bernoulli  constant  (B)  defined  as: 


B  = 


(6.6) 


In  order  to  apply  the  free  surface  boundary  conditions,  the  location  of  the  free 
surface  must  be  found,  together  with  the  potential  function,  in  each  window.  In 
order  to  locate  the  free  surface,  the  water  surface  is  defined  at  M  nodes,  located  at 
the  horizontal  location  of  the  center  of  the  array,  and  equidistantly  spaced  in  time 
throughout  the  window.  The  elevation  of  these  nodes  is  unknown,  and  will  be  sought 
as  part  of  the  solution. 


6.2  Formulation  of  the  Optimization 

The  formulation  of  the  optimization  for  the  LFI  method  as  applied  to  the  inter¬ 
pretation  of  an  array  of  pressure  measurements  has  a  great  deal  in  common  with  the 
formulation  for  an  arraj'  of  water  surface  measurements  and  the  formulation  for  a  sub¬ 
surface  pressure  record  (chapters  3  and  5).  Onh'  the  framework  is  presented,  together 
with  any  additional  information  unique  to  the  application  to  a  pressure  array. 

Observational  Equations  The  observational  equations  are  the  equations  that  con¬ 
fine  the  solution  to  fit  the  measured  records.  The  predicted  subsurface  pressure  is 
available  from  the  Bernoulli  equation,  Eq.  6.5,  which  is  applied  at  the  location  of  the 
gauges  in  the  array  to  define  the  observational  equations.  The  error  in  the  Bernoulli 
equation  is  the  difference  between  the  measured  dynamic  pressure  at  the  given  lo¬ 
cation  and  time  and  the  dynamic  pressure  predicted  from  the  kinematics  defined  by 
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the  potential  function.  Sufficient  independent  equations  are  defined  by  applying  the 
Bernoulli  equation  at  a  number  of  points  in  time  throughout  the  considered  window. 

In  order  to  specify  the  solution  to  a  system  of  equations,  there  must  be  at  least 
as  many  independent  equations  as  unknown  parameters  in  the  system.  In  order 
to  specify  the  solution  to  the  LFI  formulation  for  a  pressure  array,  the  boundary 
conditions  {f^  and  /* )  are  applied  at  each  of  the  water  surface  nodes,  yielding  2M 
equations,  and  the  Bernoulli  equation  (f^i)  is  applied  at  I  times  on  each  of  the  N 
measured  pressure  records  within  the  local  window,  yielding  N I  equations.  In  the 
single  wave  case,  there  are  4  +  J  +  M  unknown  parameters  sought  in  Eq.  6.1  (k^,  ky, 
w,  Q,  Ai . . . Aj,  and  r)i... r)M)  in  each  window,  so  that  if  (2Af  +  NI)  >  (4  +  J  +  M) 
the  solution  is  specified.  In  the  two  wave  case,  there  are  2  j  +  8  +  M  unknowns  in 
Eq.  6.2  {2kx,  2ky,  2u>,  2a,  Ai. . .  Aj,  and  rji . . .  r/M-  10  +  M  at  first  order,  14  +  M  at 
second  order,  20  +  M  at  third  order,  and  28  +  M  at  fourth  order).  The  solution  is 
specified  when  (2M  +  AT/)  >  (2Ei  +  8  +  M).  This  system  results  in  the  following 
nonlinear  least  squares  optimization. 

N  I 

minimizeO(X)  =  2/n,  tif  + 

n=l 

M  (6-7) 

y  (^5  Vci  Vmi  tmY  d"  /m  (^5  Vci 

m=l 

where: 


X  —  (kx,  ky,u,  a,  Aj, . . .  ,  Aj,  . . .  rjM) 
for  the  one  wave  case,  and: 

X  =  {^x,iiky,i,u}i,  ai,kx,2i  ky^2')^2iOc2i  Aj, . . .  ,  Aj,  %  .  tjm) 

for  the  two  wave  case.  Xn,  j/n?  a^id  Zn,p  are  the  coordinates  of  the  nth  gauge,  Xc  and  j/c 
are  the  coordinates  of  the  center  of  the  array,  and  is  the  elevation  of  water  surface 
node  m. 

As  with  the  previous  analysis,  overspecification  is  helpful  in  accommodating  the 
measurement  errors  in  field  records,  as  well  as  being  required  to  define  the  shape  of 
measured  records  and  water  surface  in  each  window. 
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6.3  Computation  Methods 

The  LFI  method  for  an  array  of  pressure  measurements  can  be  broken  down  into 
the  same  sequence  of  steps  as  with  the  analysis  of  a  water  surface  array: 

1.  Pre-processing  of  record. 

(a)  Determine  estimate  for  level  of  noise  in  the  record. 

(b)  Determine  MWL  and  subtract  hydrostatic  pressure  from  the  records. 

(c)  Determine  estimate  for  magnitude  and  direction  of  Eulerian  current. 

(d)  Define  a  set  of  continuous  records  from  each  gauge  from  the  discrete  ob¬ 
servations  by  cubic  spline  interpolation. 

(e)  Specify  spacing  of  output  locations. 

(f)  Compute  mean  zero  crossing  frequency. 

(g)  Non-dimensionalize  record  and  all  parameters. 

2.  Primary  values  of  numerical  solution  parameters  are  chosen. 

(a)  Window  width  (tq) 

(b)  Order  of  solution  (J) 

(c)  Number  of  time  samples  of  the  pressure  records  within  each  window  (7) 

(d)  Number  of  water  surface  nodes  within  each  window  (A7) 

3.  Global  solution  is  computed  on  an  entire  wave  to  provide  first  estimates  for 
local  optimization 

4.  For  each  selected  output  location,  a  window  in  the  record  is  defined,  and  an 
LFI  solution  is  computed. 

(a)  Initial  guess  for  the  optimization  determined  from  the  global  solution. 

(b)  Full  nonlinear  optimization  for  all  unknown  components  of  the  potential 
function  is  computed. 
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(c)  Results  are  checked  for  spurious  solution. 

i.  If  no  solution,  or  a  spurious  solution,  is  found,  the  solution  parameters 
are  adjusted,  and  the  optimization  repeated. 

ii.  If  a  good  solution  is  found,  progress  to  the  next  window. 

6.3.1  Pre-Processing  of  Record 

The  pre-processing  of  the  measured  records  is  essentially  the  same  as  presented  in 
the  previous  chapters.  Estimates  for  the  mean  water  level  (h)  and  Eulerian  current 
(Ux  and  Uy)  must  be  computed  to  define  the  propagation  medium.  Cubic  splines  of 
the  measured  records  are  computed  to  provide  continuous  records  for  computation, 
and  the  desired  output  locations  are  chosen.  The  records  and  all  parameters  and 
variables  are  non-dimensionalized  by  the  following  parameters:  the  mass  density  of 
water  (p),  acceleration  of  gravity  (p),  and  the  mean  zero  crossing  frequency  of  the 
measured  records  (a;^). 

6.3.2  Optimization  Procedure 

The  nonlinear  optimization  in  the  LFI  method  for  the  analysis  of  arrays  of  pressure 
measurements  is  very  similar  to  that  used  for  the  analysis  of  a  single  subsurface 
pressure  trace.  In  the  single  wave  approach,  the  solution  is  somewhat  easier.  The 
potential  function  used  by  the  single  wave  approach  is  almost  the  same  as  that  used 
with  a  point  measurement  (both  chapter  3  and  Sobey  (1992)),  with  the  addition 
of  a  directional  component  to  the  wave  number.  This  results  in  a  single  additional 
unknown,  but  the  array  of  measurements  provides  at  least  three  times  as  much  data, 
with  gauges  at  at  least  three  spatial  locations  to  specify  uniquely  the  directional 
structure  of  the  sea.  The  result  is  that  the  optimization  tends  to  converge  more 
rapidly  and  robustly  than  with  a  single  point  measurement. 

When  applying  the  method  with  the  two  wave  potential  function,  there  are  far 
more  unknowns,  and  the  optimization  once  again  becomes  somewhat  tenuous.  As 
with  the  previous  chapters,  it  is  essential  to  identify  a  good  initial  estimate  to  reduce 
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the  chances  that  the  optimization  will  converge  to  a  spurious  minimum. 


Global  Solution 

As  with  arrays  of  water  surface  measurements,  a  single  short  segment  of  the  mea¬ 
sured  pressure  records  often  may  not  contain  sufficient  information  to  identify  the 
directional  trend  of  the  wave  field.  As  a  result,  it  is  necessary  to  examine  a  larger 
segment  of  the  record,  for  example  an  entire  wave  from  crest  to  crest  or  trough  to 
trough,  to  get  a  general  sense  of  the  directional  trend  of  the  wave  field.  This  global 
solution  can  then  be  refined  to  fit  a  small  segment  very  closely.  This  approach  is  used 
to  establish  the  initial  estimate  for  the  optimization  procedure  in  each  window. 

The  procedure  for  computing  the  initial  estimate  to  the  global  solution  for  an 
array  of  pressure  measurements  is  virtually  identical  to  that  used  in  the  previous 
chapter  for  arrays  of  water  surface  measurements.  It  will  be  outlined  here.  For  more 
detail,  see  section  5.3.2. 

For  the  initial  estimate,  it  is  assumed  that  the  pressure  records  can  be  approxi¬ 
mated  with  linear  wave  theory: 


Pd  ■“  u  COS  -|-  u)  (6.8) 

for  the  single  wave  method,  and 


Pd  =  ai  cos  (A-x.iXn  -f  -  U^itm  +  Oi) 

-f  02  cos  {kr.2Tn  -f  ky^^Vn  -  +  <>2) 


(6.9) 


for  the  two  wave  method,  where: 


On  —  AnPWn 
Kn 


cosh  Kn{h-\-  Zp) 
cosh  K„h 


(6.10) 


kl  +  k] 


An  is  the  amplitude  of  the  linear  potential  function  of  the  nth  wave,  and  Zp  is  the 
elevation  of  the  pressure  array. 
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Directional  Trend  A  first  estimate  for  the  directional  trend  of  the  wave  field  is 
computed  from  the  gradients  of  the  dynamic  pressure.  A  set  of  complex  direction 
vectors  are  computed: 

+  ,6,„) 

The  spatial  gradients  are  estimated  from  finite  difference  approximation  from  the 
measured  pressures,  and  the  time  gradient  is  computed  from  a  fitted  smoothing  spline 
(de  Boor  1978)  of  the  estimated  dynamic  pressure  at  the  center  of  the  array. 

For  the  single  wave  method,  the  estimated  propagation  direction  is  the  orientation 
of  the  mean  of  the  complex  direction  vectors: 

^  =  angle  (6.12) 

For  the  two  wave  method,  the  two  dominant  propagation  directions  are  defined 
as  the  angles  of  the  means  of  the  sets  of  direction  vectors  within  tt  greater  than  and 

where:  ^ <  angle( )  <$  +  7r 

(6.13) 

where:  $  -  tt  <  angle(^m)  <  0 


less  than  the  mean  direction: 


EVequency  The  frequency  is  computed  from  the  second  derivative  of  the  pressure 
at  the  center  of  the  array. 


d^Pd,c,m 

dt^ 


(6.14) 


where  — is  computed  from  the  same  smoothed  spline  used  above.  For  the  two 
wave  method,  it  is  assumed  that  the  bimodal  sea  is  the  result  of  reflection,  so  that 
the  same  frequency  is  used  as  the  first  estimate  for  both  waves. 


Wave  Numbers  The  wave  numbers  are  estimated  from  the  linear  dispersion  rela¬ 
tion. 


(u?  —  KnUnf  =  gKn  tanh  Knh, 


^x,n  —  An  COS 


^y,n  —  An  sin  ^n  (6.15) 
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where  Un  is  the  component  of  the  Eulerian  current  in  the  direction  of  the  nth  wave. 

U.  =  y/Ul  +  Ul  COS  (tan-’  -  e\  (6.16) 


Amplitudes  and  Phases  The  amplitudes  and  phases  the  record  are  found  by 
separating  the  cosine  and  sine  components  so  that  the  equations  can  be  solved  as  a 
linear  least  squares  problem. 


?d,c  =  a  cos  [kxT  +  kyy  —  ut  +  a) 

=  6 cos  {kfX  +  kyy  —  a;/)  +  csin  {k^x  +  kyy  -  u>i) 
a  =  y/h^  -f 
o  =  arctan  (—c/6) 
for  the  single  wave  method,  and 


(6.17) 


Pd,c  =  Oi  COS  {kx^iX  +  ky  iy  -  wt  +  Oi )  -f  02  cos  {k^.^r  -f  fcj,,2y  -  wf  +  02) 
=  61  cos  {kj.^ix  -f  ky_iy  -  ut)  +  Ci  sin  (Ar^jx  +  ky,iy  -  tjjt) 

+62  cos  {kx,^^  +  ky,2y  —  u>t)  +  C2  sin  {k^.^x  -f  ky^^V  -  ‘*’0 
ai  =  y/bl  -f  cj 

02  =  •^62  +  02 

aj  =  arctan  (—Cl /61) 

Q2  =  arctan  (— C2/62) 


(6.18) 


Refining  the  Linear  Estimates  These  rough  estimates  are  refined  by  optimizing 
for  a  best  linear  wave  theory  fit  to  the  record: 

N  M 

minimize  0(X)  ^  ^  (X;  t/„,  <y„))  (6.19) 

n=l  m=l 

where: 


f  1  J/n  1  ^  COS 


X  =  (kjr.ky^a^a) 

K  =  y/ki  +  kl 
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for  the  single  wave  method,  and: 

(X;  yn^  tm)  =  Cll  COS  {kx^lX^  +  ky^ipn  “  +  CKl) 

“h  ^2  (^kx^2^n  “I"  ky^2yn  H"  ^2) 


X  {kx^i<)  ky^ij  Cki^  kx^2’)  ky^2j  ^2^)  ^2) 

K.  =  ^/K.„  +  kl, 

for  the  two  wave  method.  The  frequencies  are  determined  from  the  linear  dispersion 
relation: 


=  y/ gKn  tanh  Knh  +  KnUn 


(6.20) 


where  Un  is  defined  by  Eq.  6.16  This  optimization  results  in  a  linear  estimate  for 
the  dynamic  pressure  that  fits  the  measured  records  most  closely  in  the  segment 
considered. 

The  final  step  in  computing  the  initial  estimates  for  the  final  window-by-window 
optimization  is  to  compute  a  full  order  global  solution  to  this  large  segment  of  the 
record.  The  initial  parameters  for  this  full  order  global  optimization  are  provided  by 
the  computed  linear  fit  to  the  water  surface  records,  with  the  amplitudes  adjusted  for 
the  potential  function: 


^  _  a„  cosh  Knh 

pUJn  cosh  Kn{h  +  Zp) 

The  higher  order  Fourier  amplitudes  are  all  initially  set  to  zero. 


(6.21) 


Water  Surface  The  location  of  the  water  surface  at  M  nodes  in  the  center  of  the 
array  throughout  the  segment  is  estimated  from  the  linear  pressure  response  function 
with  stretching  (Nielsen  1989): 


Pd,  {tm)  cosh  {k  {h  -h  Pd{tm)lpg)) 
pg  cosh  k{h  -h  Zp) 


(6.22) 


where  k  is  the  mean  of  the  two  wave  numbers,  Zp  is  the  elevation  of  the  pressure 
array,  and  rjm  is  the  elevation  of  the  water  surface  node  at  the  center  of  the  array  at 
im-  Pd{tm)  is  computed  from  Eq.  6.8  or  Eq.  6.9. 
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Using  these  linear  wave  theory  estimates  as  the  first  guess,  the  full  order  optimiza¬ 
tion,  Eq.  6.7,  is  computed  to  determine  the  best  full  order  fit  to  a  global  segment  of 
the  record.  The  parameters  of  the  potential  function  computed  by  this  optimization 
are  then  used  as  the  initial  estimate  for  the  final  optimization  in  each  defined  small 
window  in  time. 

Phase  Shift  The  phase  parameters,  Oj  and  02  are  adjusted  to  accommodate  the 
change  in  the  time  reference  frame  from  the  global  to  the  local  solution; 

~  -|-  On, global  (6.23) 

where  At  is  the  difference  between  the  time  in  the  two  reference  frames. 

Nonlinear  Optimization  in  Windows 

The  established  global  solution  is  used  as  the  first  guess  for  the  parameters  in 
each  window,  and  these  parameters  are  refined  in  the  same  manner  as  the  previous 
chapters  by  solving  Eq.  6.7  with  any  standard  nonlinear  optimization  routine. 

The  resulting  solution  is  then  checked  to  see  if  is  clearly  spurious.  Spurious  solu¬ 
tions  can  be  identified  by  the  same  criteria  used  in  the  previous  chapters: 

•  Very  large  or  systematically  variable  errors. 

•  First  order  amplitude  smaller  than  one  of  the  higher  order  amplitudes. 

•  Unrealistically  large  or  small  frequency  or  wave  number. 

•  Large  discontinuity  between  the  windows  in  the  predicted  kinematics. 

If  no  solution  or  a  spurious  solution  is  found,  the  parameters  of  the  solution  are 
revised  for  another  attempt.  These  revisions  include  increasing  the  width  of  the 
window,  and  decreasing  the  order  of  the  solution.  If  neither  of  these  adjustments 
result  in  an  acceptable  solution,  the  window  is  skipped,  and  future  analysis  must  be 
interpolated  through  that  point. 

As  with  the  analysis  of  a  point  pressure  measurement,  these  adjustments  are  most 
likely  to  be  needed  in  the  long,  flat  troughs  of  shallow  water  waves,  or  near  zero 
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crossings  in  the  record.  The  additional  data  provided  by  the  array  of  measurements 
provides  for  a  more  robust  optimization  than  with  a  single  subsurface  pressure  mea¬ 
surement.  As  well  as  providing  additional  data,  the  spatial  array  provides  information 
about  the  evolution  of  the  wave  field  in  space,  helping  to  define  the  local  wave  number 
more  clearly.  This  is  in  contrast  to  the  analysis  of  a  point  measurement,  where  the 
spatial  evolution  of  the  wave  field  must  be  determined  entirely  from  the  measured 
temporal  evolution,  coupled  with  the  governing  equations.  As  a  result,  adjusting  the 
solution  parameters  is  necessary  less  often  than  with  a  single  measurement. 

6.4  Theoretical  Records 

6.4.1  Single  Wave  Records 

The  following  results  are  from  “measured”  data  generated  by  Fourier  steady  wave 
theory  (Sobey  1989),  providing  a  near-exact  solution  for  steady  waves  that  provides 
the  complete  kinematics.  The  data  used  are  a  simulation  of  data  that  might  be  col¬ 
lected  by  a  DWG-1  pressure  array  (Howell  1992).  The  DWG-1  is  a  reliable,  easy  to 
deploy  pressure  array.  The  unit  is  capable  of  including  up  to  six  independent  pressure 
transducers,  but  is  frequently  used  with  three  transducers,  arranged  in  an  1.6  m  equi¬ 
lateral  triangle  (fig.  6.1).  Three  measurements  were  chosen  for  the  theoretical  data,  as 
three  is  the  minimum  number  required  to  provide  directional  information.  Additional 
measurements  would  provide  overspecification,  and  can  easily  be  accommodated  in 
the  formulation. 

Shallow  Water  Figure  6.2  is  a  steady  shallow  wave  generated  by  18th  order  Fourier 
theory  with  the  following  parameters:  wave  height  =  3m,  period  =  10s,  water  depth 
=  5m,  direction  of  travel  10  degrees  from  the  x-axis,  and  an  opposing  Eulerian  current 
of  2m/s.  This  is  a  near  limit  wave  in  this  depth  water,  with  this  opposing  current. 
The  pressure  array  is  located  at  the  bottom.  The  parameters  of  the  LFI  solution  are: 
window  width  =  Is  (tq  =  O.IT,),  fifth  order  (J  =  5),  with  6  samples  on  the  water 
surface,  and  6  samples  on  each  of  the  pressure  records  (M  =  /  =  6),  resulting  in  30 
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wave  ray 


Figure  6.1:  Layout  of  DWG- 1  pressure  array 


equations  in  15  unknowns. 

The  solid  lines  are  the  dynamic  pressure,  water  surface,  and  kinematics  as  com¬ 
puted  from  the  Fourier  solution.  The  circles  are  the  predictions  of  the  LFI  method, 
with  each  circle  being  located  at  the  center  of  a  separate  window.  The  kinematics 
are  measured  at  a  depth  of  1.5m  below  the  surface,  just  under  the  trough.  The  LFI 
method  has  captured  the  location  of  the  water  surface  and  the  kinematics  essentially 
exactly,  including  the  pronounced  sharp  crest.  As  with  the  single  pressure  gauge 
(Sec.  3.4),  a  fairly  high  order  solution  {J  =  5)  was  necessary  to  capture  the  sharp 
crest  of  this  shallow  water  wave.  The  additional  data  provided  by  the  three  gauges 
allowed  for  a  narrower  window  than  with  the  single  measurement,  even  at  this  high 
order. 

Deep  Water  Figure  6.3  is  a  steady  deep  water  wave  generated  by  10th  order  Fourier 
theory  with  the  following  parameters:  w'av^e  height  =  10m,  period  =  10s,  water  depth 
=  100m,  direction  of  travel  10  degrees  from  the  x-axis,  and  zero  Eulerian  current.  The 
pressure  array  is  located  10m  below  the  surface.  The  parameters  of  the  LFI  solution 
are:  window  width  =  Is  (tq  =  O.lTj),  second  order(J  =  2),  with  3  samples  on  the 
water  surface,  and  3  samples  on  each  of  the  pressure  records  (M  =  /  =  3),  resulting 


Figure  6.2:  LFI  predictions  ( J  =  5)  and  exaet  kinematics  at  the  center  of  the  array 
at  z  =  — 1.5  m  for  a  steady  shallow  water  wave. 
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Figure  6.3:  LFI  predictions  (J  =  2)  and  exact  kinematics  at  the  center  of  the  array 
at  2:  =  —  5  m  for  a  steady  deep  water  wave. 


139 


in  12  equations  in  9  unknowns. 

The  LFI  method  has  captured  the  location  of  the  water  surface  and  the  kinematics 
essentially  exactly,  at  only  second  order.  Higher  order  may  be  necessary  to  define  the 
kinematics  of  a  more  irregular  record,  and  can  be  included  with  little  computational 
penalty. 

6.4.2  Records  of  Two  Intersecting  Waves 

Theoretical  pressure  records  were  generated  by  the  same  fourth  order  stokes-type 
intersecting  wave  theory  used  in  the  previous  chapter  (Ohyama,  Jeng,  and  Hsu  1995a). 
This  method  provides  the  complete  kinematics,  is  applicable  in  deep  water,  and  as¬ 
sumes  a  zero  Eulerian  current  (see  appendix  A). 

Standing  Wave  Figure  6.4  shows  the  results  of  the  method  for  a  standing  wave. 
The  parameters  of  the  Ohyama  solution  are:  period  =  10s,  first  order  amplitudes  = 
3m  (resulting  in  a  total  wave  height  of  slightly  over  6m),  propagation  directions  =  15 
and  195  degrees  from  the  x  axis,  in  deep  water.  The  pressure  array  is  located  10m 
below  the  mean  water  level.  The  LFI  method  applied  to  third  order  ( J  =  3)  finds 
the  the  water  surface  and  the  kinematics  at  the  elevation  of  3m  below  the  surface, 
just  below  the  trough  of  the  wave.  While  these  results  show  the  complete  wave,  it  is 
important  to  keep  in  mind  that  as  with  all  the  previous  results,  each  of  the  indicated 
points  is  in  the  center  of  a  separate  window,  and  was  computed  independently  of  the 
other  windows.  In  this  case,  the  standard  window  width  was  Is  (tq  =  O.IT2),  with 
a  third  order  potential  function  (J  =  3),  four  water  surface  nodes  (M  =  4)  and  six 
samples  on  the  pressure  records  (/  =  6)  distributed  equally  in  time  in  each  window, 
resulting  in  26  equations  in  24  unknowns. 

Short  Crested  Wave  Figure  6.6  shows  the  results  of  the  method  for  the  short 
crested  wave  shown  in  Figure  6.5  The  parameters  of  the  wave  are:  period  =  10s, 
first  order  amplitudes  =  3m  (resulting  in  a  wave  height  of  just  over  6m),  propagation 
directions:  30  and  150  degrees  from  the  x  axis,  in  deep  water.  The  pressure  array 


of  a  pressure  array 


kinematics  at  the  center  of  the  array 
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Figure  6.5:  Water  surface  of  short  crested  wave  at  t=0 

is  located  at  lOm  below  the  MWL.  The  LFI  method  applied  to  third  order  {J  =  3) 
again  finds  the  kinematics  for  this  wave  almost  exactly.  In  this  case,  the  standard 
window  width  was  also  Is  (tq  =  O.IT,),  with  four  water  surface  nodes  (M  =  4)  and 
six  samples  on  the  pressure  records  (7  =  6)  distributed  uniformly  in  time  in  each 
window,  resulting  in  26  equations  in  24  unknowns. 

Figure  6.7  is  the  same  short  crested  wave,  but  computed  with  the  single  wave 
method  of  the  LFI  solution.  In  this  case,  the  LFI  solution  still  matches  the  measured 
pressure  record  very  well,  as  is  virtually  always  the  case,  as  those  measurements  are 
part  of  the  optimization.  The  predicted  water  surface  is  also  fairly  accurate,  but 
with  the  predicted  crest  slightly  underestimated.  The  vertical  velocity  is  also  fairly 
accurate.  The  LFI  predictions  for  horizontal  velocities,  on  the  other  hand,  are  very 
different  from  the  Ohyama  solution.  The  resulting  discontinuities  in  the  predictions 
for  the  horizontal  velocities  make  it  clear  that  something  important  is  missing  from 
the  solution.  The  single  wave  method  is  simply  not  able  to  capture  a  short  crested 
sea  made  up  of  two  distinct  directional  components. 
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Figure  6.7:  LFI  predictions  ( J  =  3)  and  exact  kinematics  at  the  center  of  the 
at  z  =  — 3  m  for  a  short  crested  wave  computed  with  the  single  wave  method. 


144 


Choice  of  Order 

The  choice  of  order  for  the  analysis  of  arrays  of  pressure  measurements  is  essen¬ 
tially  the  same  as  for  a  single  pressure  measurement.  See  chapter  3  for  a  detailed 
discussion.  In  general,  lower  orders  are  necessary  for  deep  water  than  shallow  water, 
with  higher  than  third  order  unlikely  to  be  necessary  in  deep  water.  Shallow  water 
may  require  up  to  fifth  or  sixth  order  in  order  to  capture  near  limit  waves. 

Choice  of  Window  Width 

The  choice  of  window  width  is  also  very  similar  for  the  analysis  of  pressure  arrays 
as  it  is  for  the  analysis  of  a  single  pressure  measurement,  or  arrays  of  water  surface 
measurements.  Window  widths  of  a  minimum  of  about  1/10  the  zero  crossing  period 
are  required,  with  the  occasional  need  to  widen  the  windows  near  zero  crossings.  The 
additional  data  provided  by  the  array  of  pressure  measurements  allowed  a  solution  to 
a  very  steep  shallow  water  wave  at  high  order  with  a  window  width  of  O.lTt,  where 
a  similar  solution  required  the  window  width  to  be  doubled  to  0.2^;  when  there  was 
only  a  single  point  pressure  measurement  av'ailable  (Section  3.4). 

Number  of  Nodes  on  the  Water  Surface  and  Pressure  Records 

The  criteria  developed  in  the  previous  chapters  for  the  number  of  water  surface 
nodes  are  applicable  here.  J-f  1  points  are  required  to  specify  uniquely  a  Fourier  series 
of  order  J.  In  order  to  keep  the  order  of  the  water  surface  consistent  with  the  order 
of  the  potential  function,  at  least  M  =  J+1  water  surface  nodes  should  be  used.  The 
same  number  of  points  on  the  measured  pressure  records  is  usually  adequate.  For  the 
standing  and  short-crested  waves,  additional  equations  were  needed  to  specify  the 
solution,  so  /  =  J  -1-  3  points  on  the  pressure  records  were  used.  The  results  are  not 
highly  sensitive  to  this  parameter,  provided  sufficient  samples  are  used  to  define  the 
local  curvature  and  specify  the  system  of  equations.  Mord  points  on  the  measured 
pressure  records  will  assist  in  accommodating  noise  in  the  record. 
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6.5  Field  Measurements 

The  following  results  are  from  data  provided  by  Gary  L.  Howell  of  the  US  Army 
Corps  of  Engineers  Waterways  Experiment  Station,  Coastal  Engineering  Research 
Center  (CERC).  The  data  were  collected  as  part  of  CERC’s  Coastal  Inlets  Research 
Program.  They  were  collected  near  the  Ponce  de  Leon  inlet,  Florida,  by  a  DWG-1 
with  three  absolute  pressure  gauges,  data  sampled  at  5  Hz.  The  DWG-1  was  resting 
on  the  bottom,  so  that  the  pressure  sensor  diaphragms  were  positioned  0.21m  from 
the  bed. 

6.5.1  Field  Results 

Figure  6.8  is  a  segment  of  a  record  collected  on  August  20,  1996.  The  statistics  of 
the  record,  based  on  about  one  hour  of  record  are:  mean  wave  height  =  1.05m,  peak 
period  =  5.6s,  peak  direction  =  83°  from  true  north,  in  a  depth  of  approximately 
7.9m. 

The  top  plot  is  the  dynamic  pressure  at  the  center  of  the  array.  The  solid  line 
is  the  approximate  measured  pressure,  computed  by  linear  interpolation  from  the 
three  measured  points  in  the  array.  The  circles  are  the  values  predicted  by  the 
LFI  method.  The  other  four  plots  are  the  predicted  water  surface,  and  the  three 
orthogonal  velocities,  as  predicted  by  the  LFI  computation.  The  single  wave  method 
was  used,  with  the  following  parameters:  third  order  {J  =  3),  four  water  surface 
nodes  and  four  samples  on  the  pressure  records  (I  =  M  =  4),  and  a  window  width  of 
1.42s  (0.3T2),  resulting  in  20  equations  in  11  unknowns.  There  were  no  data  available 
about  the  Eulerian  current,  so  it  is  taken  as  zero. 

Window  Width  In  the  previous  section,  results  on  theoretically  generated  records 
indicated  the  successful  solutions  were  possible  with  quite  narrow  windows,  gener¬ 
ally  one  tenth  of  the  zero  crossing  period.  When  working  with  this  field  record,  the 
optimization  converged  with  a  window  width  that  narrow,  but  it  resulted  in  wild  fluc¬ 
tuations  in  the  predicted  propagation  direction  from  window  to  window.  Figure  6.9  is 
the  same  segment  of  the  record  as  in  figure  6.8,  computed  with  the  same  parameters. 


147 


except  that  the  window  width  used  was  much  narrower,  tq  =  O.IT^).  There  are  large 
discontinuities  in  the  horizontal  velocity  in  the  x  direction  (third  plot),  near  both 
crests.  Figure  6.10  shows  the  parameters  of  the  potential  function  from  this  solution. 
6,  where  9  =  tan  ^[ky/kx),  is  the  local  direction  of  propagation  of  the  wave.  There  is 
a  sudden  shift  in  direction  of  80°  at  the  first  crest,  and  over  100°  at  the  second  crest. 
This  would  result  in  unrealistic  predictions  of  huge  accelerations.  Figure  6.11  shows 
the  parameters  of  the  solution  computed  with  tq  —  0.3T^.  ^Vith  the  wider  window, 
the  propagation  direction  varies  smoothly  around  the  peak  direction  of  the  record, 
9  =  83°. 


Array  Size  The  need  for  a  fairly  wide  window  for  the  field  solution  is  likely  due 
to  the  small  size  of  the  pressure  gauge  array.  Compact  size  is  one  of  the  major 
advantages  of  the  DWG-1  array,  as  it  is  a  single  unit  that  is  easily  deployed,  but 
it  makes  it  susceptible  to  difficulties  with  local  analysis.  The  propagation  direction 
within  each  window  is  determined  by  the  spatial  gradients  of  the  dynamic  pressure, 
as  estimated  by  the  measurements.  In  segments  of  the  record  where  the  gradients 
are  small,  the  predicted  direction  of  propagation  can  fluctuate  a  great  deal  with  only 
small  change  in  any  of  the  measured  values.  These  changes  could  be  the  result  of 
measurement  error,  or  be  small  fluctuations  in  the  actual  dynamic  pressure.  In  either 
case,  the  large  resulting  change  in  predicted  propagation  direction  can  be  a  source  of 
error. 

The  linear  dispersion  relation  predicts  that  the  wavelength  corresponding  to  the 
peak  period  of  this  record  in  this  depth  of  water  would  be  over  42  m.  The  DWG-1 
array  is  only  1.6  m  on  a  side,  which  is  about  0.04  of  the  approximate  wavelength  of  the 
peak  period  waves.  As  a  result,  near  the  crest  of  the  wave,  where  the  water  surface 
is  close  to  horizontal,  all  three  gauges  are  near  this  point,  and  the  local  estimate 
of  the  gradient  of  the  pressure  is  very  sensitive  to  small  fluctuations.  This  effect 
is  mitigated  by  using  a  larger  window  in  time,  so  that  a  part  of  the  wave  with  a 
higher  gradient  is  part  of  the  window,  allowing  the  propagation  direction  to  be  better 
defined,  and  smoothing  out  any  measurement  errors.  As  minimum  window  size  has 
been  determined  to  be  about  1/10  of  the  zero  crossing  frequency,  an  array  size  of 
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Predicted  kinematics  at  the  center  of  the  array 
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Figure  6.9;  LFI  predictions  ( J  =  3)  of  kinematics  at  the  center  of  the  array  at  the 
water  surface  for  a  field  record,  using  a  narrow  window. 
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about  1/10  of  a  zero  crossing  wavelength  would  probably  result  in  more  stable  results 
with  the  smaller  windows. 


6.6  Discussion 

The  LFI  method  can  be  applied  to  interpretation  of  records  from  arrays  of  pressure 
measurements.  It  is  effective  for  records  from  primarily  unimodal  seas,  as  well  as  the 
bimodal  seas  that  are  generated  near  reflecting  surfaces.  The  result  of  the  analysis  is 
a  complete  description  of  the  kinematics  of  the  waves  throughout  the  water  depth,  in 
the  vicinity  of  the  array.  The  method  is  local  in  that  a  separate  potential  function  is 
used  to  describe  each  small  segment  of  the  record,  using  information  from  only  that 
segment.  This  approach  is  rational  in  that  it  seeks  to  satisfy  the  governing  equations 
of  gravity  waves,  including  the  full  nonlinear  free  surface  boundary  conditions. 

Using  a  potential  function  representing  a  single  progressive  wave,  the  method 
has  been  shown  to  be  effective  on  theoretically  generated  records  of  steady  waves  in 
both  deep  and  shallow  water.  This  formulation  is  appropriate  in  locations  far  from 
reflecting  surfaces,  where  the  sea  state  is  expected  to  be  unimodal.  When  the  method 
is  applied  with  a  potential  function  representing  two  intersecting  waves,  it  is  effective 
in  accurately  capturing  the  kinematics  of  standing  and  short  crested  waves,  as  result 
from  the  reflection  of  steady  waves  from  a  reflecting  surface,  such  as  a  sea  wall.  This 
two  wave  formulation  is  appropriate  near  such  reflecting  surfaces,  where  the  sea  state 
is  likely  to  be  bimodal. 

A  number  of  parameters  must  be  specified  to  establish  a  solution,  including  the 
order,  window  width,  number  of  water  surface  nodes  defined,  and  number  of  samples 
on  the  measured  records.  The  examples  in  this  chapter  indicate  that  the  criteria 
used  to  define  these  parameters  are  similar  to  those  developed  in  the  analysis  of  point 
pressure  records  and  arrays  of  water  surface  measurements.  Fairly  low  order  solutions 
{J  —  ^  OT  3)  are  effective  in  deep  water,  while  shallow  water  requires  higher  order 
solutions  ( J  =  5  or  6). 

Window  widths  of  1/10  the  zero  crossing  period  of  the  record  are  effective  on 
theoretical  records.  With  field  records,  however,  when  the  array  used  is  small  (less 
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than  1/10  a  typical  wavelength),  the  predictions  of  the  propagation  direction  can 
be  very  sensitive  to  subtle  fluctuations  in  the  measure  record.  In  order  to  find  a 
stable  solution,  the  window  width  must  be  increased  to  about  1/3  of  the  zero  crossing 
frequency.  This  includes  more  of  the  record  in  the  window,  allowing  the  propagation 
direction  to  be  more  clearly  defined. 

The  examples  in  this  Chapter  were  all  computed  using  MaTLAB  running  under 
the  Linux  operating  system  on  an  Intel  90MHz  Pentium  processor  based  PC.  The 
shallow  water  steady  wave  (Fig  6.2)  took  the  longest  to  compute,  about  64min.,  and 
471 X 10®  floating  point  operations  (flops).  The  steady  deep  water  wave  (Fig.  6.3)  took 
the  least  time,  12min.  and  22x10®  flops,  the  standing  wave  (Fig.  6.4)  took  16min. 
and  131x10®  flops,  the  short  crested  wave  (Fig.  6.6)  took  13min.  and  133x10®  flops, 
and  the  field  record  (Fig.  6.8)  took  17min.  and  46x10®  flops. 
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Chapter  7 
Conclusions 


This  dissertation  examines  the  problem  of  determining  the  kinematics  of  three  di¬ 
mensional  irregular  seas  from  a  variety  of  wave  measurements.  The  presented  methods 
and  results  lead  to  the  following  conclusions: 

•  Currently  used  global  methods  are  inadequate  for  determining  accurately  the 
kinematics  of  irregular  waves. 

~  Linear  superposition  does  not  satisfy  the  nonlinear  free  surface  boundary 
conditions,  resulting  in  large  errors  in  the  predicted  kinematics  near  the 
water  surface. 

—  Directional  spectral  methods  generally  disregard  the  phase  information, 
making  it  impossible  to  determine  the  detailed  kinematics  of  directional 
seas. 

“  Global  methods  of  high  enough  order  to  accurately  solve  the  free  surface 
boundary  conditions  in  three  dimensions  are  prohibitively  computationally 
intensive. 

•  The  Local  Fourier  Method  for  Irregular  waves  (LFI)  can  accurately  describe  the 
kinematics  of  two  and  three  dimensional  irregular  seas  by  satisfying  the  full  free 
surface  boundary  conditions. 
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—  The  local  nature  of  the  LFI  method  allows  it  to  satisfy  the  nonlinear  free 
surface  boundary  conditions  at  low  order  by  finding  separate  solutions  in 
each  of  a  sequence  of  small  windows  in  time. 

—  The  LFI  method  provides  an  accurate  two  dimensional  description  of  the 
kinematics  of  irregular  seas  measured  by  a  single  water  surface  gauge  or 
subsurface  pressure  gauge. 

—  The  LFI  method  provides  an  accurate  three  dimensional  description  of  the 
kinematics  of  irregular  seas  measured  by  arrays  of  water  surface  gauges  or 
arrays  of  subsurface  pressure  gauges,  including  the  bimodal  seas  resulting 
from  reflection  from  a  vertical  surface,  such  as  a  sea-wall  or  breakwater. 

—  When  working  wdth  subsurface  pressure  records,  the  LFI  method  is  limited 
in  its  capability  to  capture  the  high  frequency  components  of  the  sea  state 
that  decay  rapidly  w’ith  depth. 

—  In  order  to  accurately  capture  the  detail  of  the  me2isured  records,  the  LFI 
method  requires  very  accurate  data,  with  high  sampling  rates. 

—  The  LFI  method  can  be  adapted  to  virtually  any  arrangement  of  wave 
measuring  instruments. 

—  The  LFI  method  required  substantial,  but  not  prohibitive,  computational 
resources. 

7.1  Future  Work 

The  results  presented  indicate  the  promise  of  the  LFI  method,  and  local  methods 
in  general.  There  is  still  much  w'ork  to  be  done  before  the  method  could  be  considered 
usable  as  a  “black  box”  code. 

Much  of  this  work  revolves  around  determining  the  numerical  details  of  the  nonlin¬ 
ear  optimization.  While  the  presented  results  provide  some  guidelines  to  determining 
appropriate  parameters,  considerable  judgment  and  experimentation  with  any  given 
record  is  required.  In  the  future,  it  may  be  possible  to  establish  criteria  for  defining 
the  following  numerical  solution  parameters; 
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•  Window  width 

•  Order  of  the  potential  function 

•  Number  of  water  surface  nodes 

•  Number  of  samples  of  the  measured  record 

•  Number  of  iterations  allowed  in  the  optimization 

Some  of  the  criteria  that  might  be  used  to  help  define  these  values  are: 

•  Non-dimensional  depth 

•  Ursell  number 

•  Local  curvature  of  record 

In  principle,  the  LFI  method  may  be  extended  to  additional  arrangements  of  in¬ 
struments.  While  it  is  straightforward,  in  theory,  to  adapt  the  method  to  virtually 
any  arrangement  of  instruments,  each  different  arrangement  is  likely  to  present  a  new 
set  of  numerical  solution  difficulties.  One  of  the  sources  of  this  is  the  use  of  non¬ 
linear  optimization.  Different  optimization  problems  tend  to  be  unique,  and  require 
substantial  experience  in  order  to  determine  the  appropriate  approach. 

The  experience  gained  in  applying  the  LFI  method  to  a  variety  of  arrangements 
of  instruments  helps  demonstrate  the  need  for  appropriate  data.  Whether  this  par¬ 
ticular  method  of  analysis  is  employed,  or  any  other  nonlinear  technique,  the  need  for 
accurate  data  is  paramount.  The  design  of  a  field  data  collection  program  inherently 
includes  assumptions  as  to  how  the  resulting  data  will  be  analyzed.  Currently  used 
data  collection  programs  are  often  designed  with  the  assumption  that  the  data  will  be 
analyzed  with  the  common  methods  of  linear  spectral  analysis.  As  a  result  the  data 
may  not  be  appropriate  for  the  nonlinear  analysis  needed  to  determine  the  detailed 
kinematics  of  a  measured  wave  field. 

®f  the  considerations  that  should  be  included  in  the  design  of  a  field  mea¬ 
surement  program  are: 
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•  Array  size  appropriate  for  the  analysis  method  chosen  and  wave  field  expected. 

•  Measurement  or  prediction  of  the  local  current  field. 

•  Sampling  rate  that  captures  the  detail  of  the  record. 

•  Reduction  of  sensor  noise  and  increased  instrument  accuracy. 

As  the  methods  for  measuring  waves  and  interpreting  those  measurements  become 
more  accurate,  the  understanding  of  the  complex  processes  in  the  ocean  and  coasts 
that  are  molded  by  waves  will  increase  as  well.  This  understanding  should  lead  to  safer 
coastal  and  offshore  structures,  and  well  as  reducing  the  impact  that  these  structures 
have  on  the  sensitive  coastal  environment. 
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Appendix  A 

Intersecting  wave  theory 


A.l  Introduction 

The  theory  of  multiple  wave  interactions  presented  by  Ohyama,  Jeng,  and  Hsu 
(1995a)  has  been  used  for  two  purposes  in  this  work.  The  most  important  function 
was  to  suggest  a  form  for  a  general  three  dimensional  potential  function  representing 
interacting  waves.  The  form  chosen  is  presented  in  Chapter  4.  It  also  provided  a 
theoretical  method  for  generating  a  complete,  consistent  set  of  records  to  use  for  the 
development  and  testing  of  the  LFI  method.  This  appendix  presents  a  review  and 
detailed  critique  of  the  intersecting  wav'e  theory. 

A.2  Theoretical  Background 

Ohyama  et  al.’s  method  is  an  extension  of  well  accepted  and  verified  methods. 
It  is  essentially  a  Stokes  type  expansion  in  the  wave  steepness.  As  such,  it  is  ex¬ 
pected  have  greatest  applicability  in  deep  water.  The  flow  is  taken  to  be  irrotational 
and  incompressible,  and  thus  the  governing  equations  are  the  same  as  those  used  in 
chapter  4  and  for  most  finite  depth  irrotational  wav'e  theory: 

•  Mass  conservation  (Laplace  equation)  (Eq.  4.2) 

•  Bottom  boundary  condition  (Eq.  4.3) 
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•  Dynamic  free  surface  boundary  condition  (Eq.  4.4) 

•  Modified  kinematic  free  surface  boundary  condition  (Eq.  4.6) 

There  is  no  accommodation  for  an  Eulerian  current,  although  it  might  easily  be 
included  in  future  work.  The  non-dimensional  potential  function,  water  surface,  and 
frequencies  are  expanded  in  a  power  series  in  the  small  parameter,  (e  =  ka),  where  k 
is  a  typical  wave  number  and  a  is  a  typical  amplitude  of  the  first  order  component  of 
the  water  surface. 


(j,  =  +  0{e^)  (A.l) 

T)  =  +  e2^(2)  ^  ^3^(3)  ^4,^(4)  ^^^5^  ^^  2) 

Wi  =  cu;f )  -h  ^  +  0(e®)  (A.3) 

Where  the  77^”!  are  the  potential  function  and  water  surface  at  order  n,  and  0;^”^  is 
the  frequency  of  the  ith  wave  at  order  n.  The  steepness  of  all  the  intersecting  waves 
is  taken  to  be  of  the  same  order.  Modern  Stokes  theories  have  found  it  necessary 
to  use  an  expansion  parameter  based  on  the  physical  wave  height,  rather  than  the 
amplitude  of  the  first  order  component  (Fenton  1990).  This  is  not  feasible  with 
intersecting  waves,  as  there  is  no  clearly  defined  wave  height.  While  not  ideal,  Stokes 
methods  based  on  this  first  order  expansion  parameter  can  be  effective  (Skjelbreia 
and  Hendrickson  1962),  as  long  as  they  are  algebraically  correct  and  not  pushed  to 
near  limit  waves  or  shallow  water. 

Care  must  be  taken  when  computing  a  given  order  solution  when  multiple  parts 
of  the  solution  are  all  expanded  in  a  perturbation  series.  In  this  case,  the  potential 
function,  water  surface,  and  frequencies  are  expanded  separately,  but  the  potential 
function  and  water  surface  are  dependent  on  the  full  order  frequency.  Ideally,  a 
given  order  component  will  depend  only  upon  parameters  of  the  same  order.  In  this 
case,  however,  the  full  order  frequency  must  first  be  computed,  and  then  used  as 
the  frequency  at  all  orders.  If  the  order  of  the  frequency  used  was  matched  to  the 
order  of  77  or  0  being  computed,  the  resulting  solutions  would  be  out  of  phase,  with. 
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for  example,  having  a  different  frequency  to  It  could  be  argued  that  this 
approach  would  be  accurate  to  second  order,  but  the  components  would  be  shifted 
out  of  phase  as  time  progressed.  It  is  more  accurate  to  compute  the  results  using  the 
full  order  frequencies  at  all  orders  of  <f>  and  rj. 

The  first  order  solution  for  N  intersecting  waves  is  the  familiar  superposition  of 
linear  waves: 


N 


«=1 


a,  cosh  k{h  +  z) 
oji  cosh  kh 


sin  {kj-T  +  kyy  —  ojt  +  q) 


(A.4) 


9 


(1) 


A’ 


a,  cos  (kj.  +  kyy  -  wt  +  o) 


(A.5) 


=  \/gk,  tanh  /r./j,  (A.6) 

The  higher  order  solutions  are  computed  by  applying  the  Laplace  equation,  the 
bottom  boundary  condition,  and  the  free  surface  boundary  conditions,  expanded  in 
Taylor  series  about  the  mean  water  level.  The  algebra  and  the  resulting  solutions 
are  long  and  extremely  tedious,  are  presented  in  Ohyama  et  al.  (1995a),  and  they 
will  not  be  repeated  here.  It  should  be  noted  that  the  frequency  was  only  solved  to 
third  order,  as  the  fourth  order  frequency  would  require  a  solution  to  the  fifth  order 
potential  function,  which  w'as  not  presented. 

A. 3  Analytical  Verification 

Ohyama  et  al.  provided  a  number  of  comparisons  with  other  wave  theories.  When 
the  method  is  applied  for  a  single  wave  component,  it  represents  a  steady  wave,  and  as 
such  should  be  expected  to  provide  the  same  solution  as  conventional  Stokes  theory. 
The  authors  presented  a  comparison  with  the  Fenton  (1985a)  solution.  The  Fenton 
solution  is  based  on  an  expansion  parameter  based  on  the  wav^e  height,  rather  than 
the  first  order  amplitude:  c  =  kH/2.  By  soKdng  for  c  in  terms  of  the  c.  The  two 
solutions  were  found  to  coincide  exactly. 
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When  applied  for  two  wave  components  of  the  same  wavelength,  same  amplitude, 
and  opposing  directions,  the  solution  is  a  standing  wave.  The  resulting  solution 
coincides  almost  exactly  with  a  standing  wave  solution  given  by  Chen  (1988),  except 
for  one  coefficient.  The  authors  suggest  that  the  slight  discrepancy  is  the  result  of  a 
typographical  error  in  Chen’s  paper. 

Ohyama  et  al.  also  indicated  that  they  compared  their  solution  with  those  for 
short-crested  seas  (Hsu  1990;  Hsu  and  Chen  1992;  Fenton  1985c)  and  once  again 
found  essentially  perfect  agreement. 

A. 4  Numerical  Verification 

The  above  theoretical  verifications  indicate  that  the  method  is  correct  when  re¬ 
duced  to  two  particular  simplifications.  It  remains  to  verify  that  the  method  is 
accurate  when  applied  to  a  more  complex  situation,  with  a  number  of  intersecting 
waves  of  different  heights,  directions  and  wavelengths. 

The  actual  computation  of  the  general  solution  is  quite  complicated  and  contains 
a  very  large  number  of  coefficients  that  must  be  computed.  Despite  this  complexity, 
once  the  code  is  written  and  debugged,  it  is  a  trivial  matter  for  modern  computers 
to  calculate  the  full  solution.  It  is  not,  however,  a  trivial  matter  to  write  the  code 
without  errors  or,  indeed,  to  find  the  likely  typographical  errors  in  the  published 
paper.  A  number  of  these  typographical  errors  have  been  presented  in  the  published 
erratum  (Ohyama  et  al.  1995b). 

An  additional  error  has  been  found  and  confirmed  by  personal  communication 
with  the  authors,  pp.  55  Eq.  47  in  Ohyama,  Jeng,  and  Hsu  (1995a)  should  read: 

C  =  I  £^^220  +  £^/5420  I  COS  2kx  -|-  £*0440  COS  4kx  -I-  1 2e  -b  £^031 1 1  cos  kx  cos  ut 
+  ^^^0222  +  £**^422 1  cos  2kx  COS  2ujt  -f  £^0313  COS  kx  COS  Sujf 
+£^^331  cos  Zkx  cos  ut  +  £^0333  COS  Zkx  cos  Zu>t  +  £^^0424  COS  2kx  cos  4u;t 
+£*0442  cos  4kx  cos  2u}t  +  £^0444  COS  4kx  cos  4ujt  -b  0{£^} 

As  a  result  of  the  complexity  of  the  computational  code,  it  is  advisable  to  use  numer- 
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ical  methods  to  verify  the  validity  of  the  solution  and  resulting  computational  code. 
The  following  results  were  produced  by  a  Fortran  code  derived  from  the  code  written 
and  generously  provided  by  Takumi  Ohyama  that  was  used  to  compute  the  water 
surface  plots  given  in  (Ohyama  et  al.  1995a).  The  code  was  expanded  to  compute 
the  full  kinematics. 

A.4.1  Richardson  Extrapolation 

Fenton  (1985a)  presented  a  variation  of  Richardson  extrapolation  that  can  be  used 
for  numerical  verification  of  wave  theory  code.  Richardson  extrapolation  (also  known 
as  extrapolation  to  the  limit)  is  a  method  traditionally  used  to  increase  the  accuracy 
of  finite  difference  computations  (Salvadori  and  Baron  1961;  Roach  1976).  In  Fenton’s 
adaptation,  the  method  is  used  to  simultaneously  check  all  the  terms  in  the  solution, 
and  to  provide  an  estimate  of  the  order  of  the  accuracy  of  the  result. 

The  form  of  Ohyama  et  al.’s  solution  solves  the  field  equation  (Laplace)  and  the 
bottom  boundary  condition  exactly.  It  is  expected  to  satisfy  the  free  surface  boundary 
conditions  to  the  order  it  is  applied.  The  following  method  calculates  the  order  of 
accuracy  of  the  free  surface  boundary  conditions.  Classical  Richardson  extrapolation 
is  used  to  evaluate  the  error  of  a  solution  at  a  given  point  in  time  and  space.  By 
taking  advantage  of  the  periodicity  of  wave  motion,  Fenton’s  variation  can  be  used 
to  evaluate  the  solution  throughout  a  complete  wave. 

Single  Wave 

The  solution  for  a  single  wave  is  periodic  and  symmetric  about  the  crest,  and  thus 
the  error  in  each  of  the  boundary  conditions  can  be  expanded  in  a  single  sided  Fourier 
series: 

£(e)  =  If 

k=0 


(A.7) 
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Fourier  Component,  (k) 

0 

1 

2 

3 

4 

5 

6 

Order  of  KFSBC  error 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 

Order  of  DFSBC  error 

3.0 

2.2 

2.0 

2.0 

2.0 

2.0 

2.0 

Table  A.l:  Order  of  errors  for  a  single  wave  computed  to  order  1 


Fourier  Component,  (k) 

0 

1 

2 

3 

4 

5 

6 

Order  of  KFSBC  error 

3.0 

3.0 

3.0 

3.0 

3.0 

3.0 

3.0 

Order  of  DFSBC  error 

2.9 

3.0 

3.0 

3.0 

3.0 

3.0 

3.0 

Table  A.2:  Order  of  errors  for  a  single  wave  computed  to  order  2 


The  magnitude  of  this  error,  and  thus  each  of  the  coefficients  can  be  expanded  in  a 
first  order  Taylor  series. 


C'/bCeO  =  Mk  +  o  (el)  (A.8) 

These  coefficients  are  a  function  of  the  small  parameter,  c  and  thus  are  a  function  of 
c”  at  nth  order, 


=  +  (A.9) 

When  the  discreet  transforms  of  the  errors  are  computed  for  two  different  values  of  e, 
an  approximation  for  the  njt  can  be  computed,  giving  an  approximation  to  the  order 
of  the  errors  of  each  harmonic  in  each  boundary  condition. 

„  _log(C4e2)/C.(€x))  .  .  . 

log(£2/€0  (A.IO) 

Table  A.l  presents  the  results  for  a  single  wave,  computed  to  order  1  in  deep  water 
{kh  —  100,  ei  =  0.01,  €3  =  0.02).  These  results  show  that  the  error  is  of  at  least 
second  order  in  both  of  the  boundary  conditions,  indicating  that  the  computation  is 
accurate  to  first  order. 

Tables  A.2  through  A.4  present  the  results  for  the  same  wave,  computed  to  sec¬ 
ond  through  fourth  order.  These  results  indicate  that  the  error  in  both  boundary 
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Fourier  Component,  (k) 

0 

1 

2 

3 

4 

5 

6 

Order  of  KFSBC  error 

4.0 

4.0 

4.0 

4.0 

4.0 

4.0 

4.0 

Order  of  DFSBC  error 

4.0 

3.9 

4.0 

4.1 

4.0 

4.0 

4.0 

Table  A.3:  Order  of  errors  for  a  single  wave  computed  to  order  3 


Fourier  Component,  (k) 

0 

1 

2 

3 

4 

5 

6 

Order  of  KFSBC  error 

5.0 

5.0 

5.0 

5.0 

5.0 

5.0 

5.0 

Order  of  DFSBC  error 

4.8 

5.0 

5.0 

5.0 

5.0 

5.0 

5.0 

Table  A.4:  Order  of  errors  for  a  single  wave  computed  to  order  4 


conditions  is  of  at  least  one  order  higher  than  the  order  used  to  compute  the  results. 
These  results  serve  to  confirm  the  efficacy  of  the  Richardson  extrapolation  method, 
as  well  as  the  accuracy  of  the  w'ave  theory  up  to  fourth  order. 

Standing  Wave 

Ohyama  et  al.’s  method  can  also  be  used  to  compute  a  standing  wave  by  defining 
two  waves  of  the  same  w'avelength  and  amplitude,  and  opposing  directions.  In  this 
case,  the  wave  is  not  steady,  and  the  errors  in  the  free  surface  boundary  conditions 
must  be  considered  for  a  full  period  in  time  and  full  wavelength  in  space.  The  ex¬ 
trapolation  to  the  limit  method  can  be  extended  to  this  ca.se  by  expanding  the  free 
surface  errors  in  a  two  dimensional  Fourier  series  (Newland  1993). 

oo  oo 

=  E  E  (A.U) 

it=0  /=0 

An  estimate  of  the  order  of  the  errors  can  then  be  computed  in  a  similar  manner  as  Eq. 
A. 10.  The  results  of  this  computation  for  a  standing  wave  in  deep  water  {kh  =  100, 
Cl  =  0.01,  C2  =  0.02)  are  given  in  Tables  A. 5  and  A. 6.  There  are  no  values  given 
for  the  zeroth  order  in  time  of  the  kinematic  boundary  condition  because  C(c)it,o  for 
both  values  of  c  are  zero  to  the  precision  of  the  calculation,  so  the  values  computed 
are  spurious.  These  results  indicate  that  Ohyama  et  al.’s  method  is  also  accurate  to 
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Fourier  Component,  (k) 

0 

1 

2 

3 

4 

5 

6 

1=0 

- 

- 

- 

- 

- 

- 

- 

1=1 

5.0 

5.0 

5.0 

5.0 

5.0 

5.0 

5.0 

1=2 

5.5 

5.5 

5.8 

5.1 

6.2 

5.6 

6.2 

1=3 

5.0 

5.0 

5.0 

5.0 

5.0 

5.0 

5.1 

1=4 

4.6 

4.9 

6.1 

5.0 

5.8 

5.2 

5.3 

1=5 

4.9 

5.1 

4.8 

5.0 

5.1 

5.0 

5.0 

1=6 

5.0 

5.0 

5.1 

5.0 

5.0 

5.0 

5.3 

Table  A.5:  Order  of  KFSBC  errors  for  a  standing  wave  computed  to  order  4 


Fourier  Component,  (k) 

0 

1 

2 

3 

4 

5 

6 

1=0 

3.6 

5.0 

5.5 

5.0 

- 

5.0 

5.1 

1=1 

3.4 

5.0 

4.1 

5.0 

4.9 

5.0 

5.0 

1=2 

5.6 

5.0 

5.7 

5.0 

6.7 

5.4 

5.3 

1=3 

5.1 

5.0 

5.1 

5.0 

5.0 

5.0 

5.0 

1=4 

6.1 

5.2 

6.3 

5.0 

5.8 

5.5 

5.3 

1=5 

4.9 

5.0 

5.2 

5.1 

5.2 

5.0 

4.9 

1=6 

5.3 

5.0 

5.9 

5.0 

4.7 

5.0 

5.3 

Table  A.6:  Order  of  DFSBC  errors  for  a  standing  wave  computed  to  order  4 


approximately  fifth  order  for  a  standing  wave. 

A.4.2  More  Complex  Seas 

Fenton’s  method  relies  on  the  periodicity  of  the  solution  to  be  tested.  If  the 
solution  is  not  periodic  in  space  and  time,  then  the  errors  in  the  boundary  conditions 
will  not  be  periodic,  and  the  Fourier  expansion  cannot  be  strictly  applied.  Fourier 
coefficients  could  still  be  computed  by  assuming  periodicity,  but  the  discontinuity 
generated  by  the  assumption  of  periodicity  would  generate  spurious  results.  These 
could  be  minimized  by  taking  a  large  number  of  points  over  a  large  range  in  time  and 
space,  but  this  would  become  very  computationally  intensive. 

The  accuracy  of  any  approximate  solution  can  be  measured  by  examining  the 


r 


174 


order  (n) 

c” 

KFSBC  error 

DFSBC  error 

1 

T 

o 

X 

1.92  X  10-" 

1.03  X  10-" 

2 

1  X  10-'' 

3.43  X  lO-'* 

1.14  X  10-3 

3 

1  X  10"'* 

1.29  X  10-3 

9.33  X  10-“ 

4 

1  X  10-“ 

3.42  X  10-“ 

1.40  X  10-“ 

Table  A. 7;  RMS  errors  in  free  surface  boundary  conditions  for  a  short-crested  sea 


errors  in  the  governing  equations  resulting  from  the  solution.  Fenton’s  method  is  a 
particularly  powerful  method  for  examining  these  errors,  in  that  it  considers  an  entire 
wavelength  and  period  at  once  and  gives  an  estimation  of  the  order  of  accuracy  of  the 
solution.  A  simpler  approach  can  also  be  useful  and  be  universally  applied.  In  the 
case  of  Ohyama  et  al.’s  method,  the  field  equation  and  bottom  boundary  condition 
are  exactly  satisfied  at  all  orders  of  approximation.  If  the  solution  is  accurate,  the 
errors  in  the  free  surface  boundary  conditions  will  decrease  with  each  increase  in  the 
order  of  approximation.  While  this  method  does  not  guarantee  that  the  solution  is 
correct,  most  errors  in  the  algebra  or  the  computational  code  will  result  in  a  decrease 
in  the  accuracy  of  the  result  at  the  order  in  which  the  error  occurs. 

Table  A. 7  contains  the  errors  in  the  free  surface  boundary  conditions  at  orders 
one  through  four  for  a  short  crested  sea.  This  sea  is  created  by  two  waves  of  the  same 
height  and  wavelength  (c  =  0.1,  kh  =  100)  intersecting  at  an  angle  of  30  degrees. 
The  values  given  are  the  root  mean  square  of  the  errors  computed  by  averaging  over 
a  grid  of  points  one  w'avelength  in  both  horizontal  directions,  and  over  one  period.  In 
this  case  the  errors  in  the  free  surface  boundary  conditions  decrease  with  increasing 
order  of  solution.  It  should  be  noted,  however,  that  the  errors  do  not  decrease  as 
much  as  c".  This  might  be  partially  due  to  the  fact  that  the  frequency  is  expressed  to 
only  second  order.  Including  the  fourth  order  frequency  would  probably  reduce  the 
error  at  fourth  order.  This  behavior  warrants  closer  examination,  particularly  if  the 
method  were  to  be  expanded  to  higher  order. 
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A. 5  Depths  of  Validity 

As  the  method  is  very  similar  to  Stokes  theory,  it  is  expected  to  be  valid  in  a  similar 
range  of  depths.  Stokes  theory  is  valid  if  both  e  and  the  parameter,  tl{khf,  similar 
to  the  Ursell  number,  are  small  (Fenton  1985a).  This  indicates  that  the  method  has 
optimal  applicability  in  deep  water,  and  should  not  be  applied  in  shallow  water. 
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13.  (CpndudMl). 


The  LFI  method  was  introduced  by  Sobey  (1992)  as  a  two  dimensional  method  for  interpreting  a  single  water 
surface  measurement  of  inegular  waves.  In  diis  dissertation,  the  two  dimensional  mediod  is  expanded  to  include  the 
analysis  of  subsurface  pressure  records.  The  method  is  dten  extended  to  the  prediction  of  directional  irregular  wave 
kinematics  from  the  measurements  of  arrays  of  instruments,  including  wave  staffs  and  subsurface  pressure  gauges. 

In  diis  extended  analysis,  the  LFI  method  includes  nonlmear  wave-wave  interactions,  in  addition  to  die  self-wave 
interactions  of  the  two  dimensional  form.  Results  of  the  analysis  (« theoretical,  laboratory  and  field  records  are 
provided. 


